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Abstract 


Branch points arise in optical transmissions due to strong atmospheric turbulence, long 
propagation paths, or a combination of both. Unfortunately, these conditions are very often 
present in desired operational scenarios for laser weapon systems, optical communication, 
and covert imaging, which suffer greatly when traditional adaptive optics systems either 
cannot sense branch points or implement non-optimal methods for sensing and correcting 
branch points. Previous research by Pellizzari presented a thorough analysis of various 
novel branch point tolerant reconstructors in the absence of noise. In this research a realistic 
model of the Air Force Institute of Technology’s adaptive optics system is developed and 
used for the first realistic tests of these branch point tolerant reconstructors. Utilizing 
a self-referencing interferometer as the high-order wavefront sensor-this type of sensor 
being theoretically immune to scintillation-this effort has extended previous theoretical 
work by adding realistic noise effects to the SRI’s measurements before reconstructing the 
wavefronts and applying the control law. Not only is the adaptive optics correction shown 
to outperform the traditional techniques by as much as 126%, but several new theories 
and refinements to existing theories were discovered along the way. This study provides 
a foundation to guide hardware implementation in the future, where noise effects will be 
present. 
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IMPLEMENTATION OF BRANCH-POINT-TOLERANT WAVEFRONT 


RECONSTRUCTOR FOR STRONG TURBULENCE COMPENSATION 


I. Introduction 

W ithin the Department of Defense (DoD), directed energy (DE) capabilities 
have been identified as one of “the key ‘game-changing’ technology-enabled 
capabilities that enters service during [the next two decades]” [44]. DE spans a wide 
range of potential applications yet has many technological challenges still to be overcome. 
Applications of particular interest to the DoD are space situational awareness, laser 
communication, and laser weapon systems [44], Currently, a major obstacle in each 
of these applications is accounting for the optical aberrations induced by the random 
nature of Earth’s atmosphere. For a point of reference, these same atmospheric effects 
are what prompted NASA to spend billions of dollars to design and maintain the Hubble 
Space Telescope. To deal with these atmospheric aberrations, there are only a handful of 
possibilities. Adaptive optics (AO) utilizes a wavefront sensor (WFS) and reconstructor to 
estimate the optical distortions on some reference beam traveling through the atmospheric 
path of interest. Subsequently, a wavefront corrector, such as a deformable mirror 
(DM), removes the estimated distortions. While there are other methods for dealing with 
Earth’s random effects on light, AO is unique in providing real-time, improved imaging 
performance, which is essential for the above applications of interest. 

Current AO techniques only provide good correction capability for limited operational 
regimes. As the optical propagation paths lengthen and/or the atmospheric turbulence 
becomes stronger, the propagated light begins to constructively and destructively interfere 
with itself [5]. In these circumstances, traditional wavefront reconstruction techniques 
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used to drive the optical correction process are unable to completely estimate the optical 
aberrations. These situations are of definite operational interest to the DoD, but most 
current AO techniques lose effectiveness in these more challenging situations [31,39,42]. 
The research presented here is aimed at extending the operational regime of AO to 
offer improved correction capability when optical interference begins to severely degrade 
performance of traditional AO techniques. 

1.1 Problem Statement 

Branch points arise in optical transmissions due to strong atmospheric turbulence, 
long propagation paths, or a combination of both. Unfortunately, these conditions are 
very often present in desired operational scenarios for laser weapon systems, optical 
communication, and covert imaging, which suffer greatly when traditional adaptive 
optics systems either cannot sense branch points or implement non-optimal methods for 
estimating and correcting branch points. The overall goal for this research is to transition 
previous theoretical work in optimized wavefront reconstruction for strong turbulence 
AO to an initial implementation in hardware. To bridge the transition, a high-fidelity 
computer model is developed and used in testing the new algorithms prior to hardware 
implementation. The work in simulation is shown to be invaluable in further investigating 
the behavior of the new wavefront reconstruction algorithms with many unanticipated 
lessons and discoveries. 

1.2 Key Results 

This research primarily investigated the behavior and properties of post-processing 
congruence operation (PCO) based wavefront reconstructors. The following list summa¬ 
rizes the contributions made in this research towards bringing the PCO-based reconstructor 
to an implementable state. 
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• Demonstrated greatly improved strong turbulence AO performance in terms of both 
mean Strehl and Strehl variance. The AO performance gains using a PCO wavefront 
reconstructor over the traditional least squares reconstructor in terms of mean field 
estimated Strehl ratio and Strehl ratio normali z ed standard deviation are summarized 
in Table 1.1. Mean Strehl is a metric of absolute AO correction performance (higher 
is better), and the Strehl standard deviation is a metric of correction stability (lower 
is better). These results were presented at the 2012 IEEE Aerospace conference. 


Table 1.1: Aggregate AO performance gains of all PCO-based wavefront reconstructors 
over the traditional least squares reconstructor. 



Mean Strehl 

Normalized Strehl Standard Deviation 

0.04 

0.2% 

-2.4% 

0.5 

19.1% 

-50.6% 

1 

126.4% 

-57.6% 


• Identified cause of fades from PCO-based AO and demonstrated a potential solution. 
When in closed-loop operation in strong turbulence, PCO algorithms can cause 
short-term fades as h opt varies. These drops in Strehl are not caused by the higher- 
order control law as much as different tilt components associated with each PCO 
rotational realization. The tilt is removed from the higher-order corrections and 
instantly dumped on the tracking system; causing the fades. One solution was found 
by allowing tilt in the higher-order corrections with conservative controller gains. 
This prevented the fades with a slight tradeoff in mean Strehl ratio. Unfortunately, a 
realistic system with slight misalignments would require a zero-DC-gain control law 
to implement a similar solution. Additionally, this puts higher stroke requirements 
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on the DM when tilt is being included. With this new understanding future work in 
control design could help mitigate these fades based on the application requirements. 

Demonstrated and developed theory on improved PCO performance in weak 
turbulence, low SNR regimes. In weak turbulence, where no branch points are 
present, a consistent performance improvement of the PCO reconstructors over 
traditional least squares was found. Although the difference is small, this was 
completely unexpected based on previous concerns of noise corrupting the PCO 
process. A new theory proposed by the author shows in high-noise conditions, 
least squares reconstructors can introduce phase anomalies to the AO corrections. 
However, the PCO process removes these phase discrepancies introduced by the least 
squares reconstruction. These findings will be presented at the 2012 SPIE Optics and 
Photonics conference. 

Developed a set of theories regarding phase anomalies in high wavefront gradient 
conditions. When localized wavefront phases being reconstructed from an SRI are 
greater than n radians between adjacent subapertures, the least squares reconstructor 
will alias, leading to a phase anomaly in the reconstructed field. This can be mitigated 
using a higher density of SRI subapertures, at the cost of SNR. Alternatively, initial 
results show the PCO corrects the phase anomalies introduced by the least squares 
reconstruction. 

Formulated an optimized data flow for an AO system using PCO. Started by 
diagnosing coupling between the high-order and tracking control systems when using 
a PCO-based wavefront reconstructor. Since PCO reconstruction is not a single 
matrix operation, decoupling tilt from the higher-order commands required a new 
data flow process. The resulting process, illustrated in Fig. 3.4, minimizes matrix 
multiplies while maximizing AO effectiveness. 
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• Conducted first beam projection simulations using PCO wavefront reconstruction 
techniques. Conducted at the end of this research, this establishes a foundation 
for future efforts to investigate these novel reconstruction methods with respect 
to beam projection applications. Although the initial simulation results identified 
needed changes, initial observations indicate promise over traditional least squares 
reconstruction. The identified flaws have since been fixed and follow-up simulation 
results will be presented at the 2012 DEPS Beam Control Conference. 

• Developed an extensible AFIT AO system model for future research. The model 
is configurable to varying levels of realism and AO capabilities. Allows isolation 
of the research problem at hand without the difficulties of hardware work and the 
assumptions of analytical work. This includes the first modeling of realistic fiber 
coupling to our specific SRI optical setup - crucial to realistic noise modeling. 

1.3 Thesis Overview 

Chapter II provides an overview of traditional AO and the optical effects as a result 
of the randomized atmosphere. This is followed by an introduction to branch points and 
the previous theoretical work in developing an optimized branch-point tolerant wavefront 
reconstructor. Chapter III describes a high-fidelity computer simulation developed to 
model AFIT’s AO system. This model exposed new insights, crucial to a successful 
implementation of the new reconstruction techniques. These insights, discussed in Ch. IV, 
offer explanations to a number of previously posed questions, and a new extension of the 
theory is presented for beam projection applications. Finally, Ch. V presents a summary of 
the challenges overcome, key contributions, and suggested future work from this research 
effort. 
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II. Background and Related Research 


This chapter introduces the current state of AO and previous theoretical methods 
for extending AO operability into regimes with scintillation present. An overview of 
a typical AO system is given, along with two WFS types - each with unique strengths 
and weaknesses. Following this, the detrimental effects of scintillation are explained 
with respect to AO. This chapter concludes by discussing previous research in specialized 
wavefront reconstruction techniques that attempt to improve upon traditional reconstructors 
when branch points are present in the received beacon wavefronts. 

2.1 Traditional Adaptive Optics 

As light propagates through Earths atmosphere, stochastic variations in temperature 
and pressure cause the light to both diffract and refract [5], This causes beam wander, beam 
spreading, and phase fluctuations. A simple example of this can be observed in the apparent 
twinkle of stars at night. The Hubble Space Telescope was built to operate beyond the 
random atmosphere, where the stars appear as sharp objects. For many applications, these 
atmospheric effects severely limit performance. For instance, beam projection systems 
rely on efficient delivery of light to a small target area, and optical communication must 
maintain a reliable optical link. One remedy is to estimate the atmospheric distortions 
and correct for them in real time. This approach is called AO, utilizing a WFS and 
reconstructor to estimate the optical distortions on some reference beam followed by a 
wavefront corrector to mitigate the estimated distortions. 

2.1.1 Atmospheric Turbulence. 

Within the atmosphere, light is distorted as it encounters spatially and temporally 
randomized variations in the air’s index of refraction. The statistical behavior of these 
random variations were first described in the 1940’s by Kolmogorov, when he used physical 


6 




Figure 2.1: Turbulent eddies shown of size L 0 down to /q. 


insights to describe the velocity of air within the atmosphere [5]. Energy cascade theory 
says that kinetic energy injected to a viscous fluid creates large sized turbulent eddies. 
These large eddies break down into progressively smaller eddies before the energy is 
eventually transferred into heat. The dimensions of these eddies are described by the 
inertial subrange, where eddie sizes fall between the inner scale, / 0 , and outer scale, L 0 . 
Within the inertial subrange, the turbulence can be considered statistically homogenous 
and isotropic [5]. Figure 2.1 provides a visual depiction of the energy cascade theory as 
kinetic energy is converted down into thermal energy. 
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With this physical foundation, Kolmogorov was able to formulate a one dimensional 
structure function for wind velocities within the inertial subrange [5], From this, one 
can derive three dimensional power spectrum of the velocity fluctuations within the 
atmosphere [17]: 

<S>rr(k) = 0.066e 2/ V 11/3 , (2.1) 

where k is the scalar spatial frequency (in units of rad/m) and e is the average energy 
dissipation rate of the atmosphere. By defining C\ as the velocity structure constant, 

C 2 = 2e 2/3 , (2.2) 

Eq. (2.1) can be simplified to 

<S>sr(k) = 0.033C 2 *- 11/3 . (2.3) 

The velocity structure constant, with units of m 4/3 /s 2 , provides a quantitative evaluation of 
the energy within the atmospheric turbulence and is actually based in the formulation of 
the velocity structure function [5]. 

From Kolmogorov’s work, Obukhov [28] and separately Corrsin [8] were able 
to extended his statistical model of atmospheric velocity variations, first to potential 
temperature, and later to index of refraction variations. This important final leap directly 
relates the phase aberrations accumulated by light to the randomized atmosphere. The 
resulting three dimensional refractive index power spectrum, known as the Kolmogorov 
spectrum, is given by 

® n (K) = 0.033C 2 k“ 11/3 , (2.4) 

where n is the index of refraction and C 2 is the index of refraction structure constant. C 2 is 
a physical quantity that can be measured for a given location, time of day, and altitude. 
Several models of C 2 exist that can be used for simulations of a desired scenario [5]. 
Figure 2.2 shows the commonly used Hufnagel-Valley model, along with many of the 
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Figure 2.2: The Hufnagel-Valley C\ profile model. Figure from [23]. 

physical atmospheric layers that drive its behavior. For this research, a horizontal path 
was assumed - implying a constant C\ along the propagation path. 

It is important to note the power spectral models for velocity, temperature, and 
index of refraction variations all follow similar behavior - easily seen by comparing 
Eqs. (2.3) and (2.4). However, as presented thus far, effects outside the inertial subrange 
have been ignored. The Kolmogorov spectrum actually assumes an infinitely large outer- 
scale and infinitely small inner-scale. This is an assumption made for mathematical 
simplicity, but does not match the true atmospheric fluidics. Hill developed a numerical 
model based on a more detailed hydrodynamic analysis [20], later adapted to analytic form 
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by Andrews [4], given by 


[ / k\ Ik \ 7/6 l exp (-K 2 /kt) 

1 + 1.802 - -0.254 - ---(2.5) 

Vi) Vi) J(*2 + ^) 11/6 

0 < k < oo, Kj = 3.3// 0 , K 0 = 2n/L 0 . (2.6) 

This is sometimes called the modified atmospheric spectrum, where the math to the right 
of C\ has been added to Eq. (2.4) to match the results found by Hill. Most notably, this 
spectrum model includes a spectral bump at wavenumbers near the inner-scale size. As 
the inner-scale size decreases, light begins to diffract at the smaller eddies [5]. With the 
Hill spectrum’s increased content of small eddie sizes, scintillation effects are seen earlier 
than a similar propagation using a Kolmogorov spectrum. For this research, the modified 
spectrum was used in all numerically propagated optical fields. 

Mathematically, Maxwell’s Equations govern all electromagnetic interactions. In 
optics, these equations are typically simplified by making three key assumptions: 

1. Negligible backscattering and polarization effects. 

2. The Markov approximation. Assumes the atmospheric refractive index variations are 
delta correlated in the direction of propagation. 

3. The paraxial approximation. Assumes all light rays of interest travel at angles very 
near the optical axis. 

From these simplifications, several methods exist for applying Maxwell’s equations to 
propagating an initial field through a random medium. For regimes of weak turbulence, 
one such method is the Rytov approximation [5]: 

U( r, L) = C/ 0 (r, z = L) exp [Y(r, L )], (2.7) 

where the initial field is propagated a distance L, Uo(r, L ) is the field propagated through 
vacuum, and it is modified by the complex phase perturbation, 'P(r, L). Further, the 
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perturbation can be split into amplitude and phase components [38] 


x ¥= X + i<t>. 


(2.8) 


(2.9) 


Several statistical parameters of the propagated field U(r,L ) can be found using the 
Rytov approximation described above. First, define the field’s mutual coherence function 
(MCF), 

T (r, r',L) = <f/(r)f/*(r')> 

= [/(r)[/*(r , )(exp [^(r)^*(r')]>. 

By taking the modulus of the MCF, the field’s spatial coherence can be found - also called 
the degree of coherence: 

F (r, r',L)|_ 

(2.10) 


r(r,r \L) = 


Vr (r, r, L) T (r\ r', L) 


1 


where D(r,r',L ) is the wave structure function (WSF) [17]. Assuming homogeneity, 
isotropy, and a spherically expanding source, the WSF is 

D sph (p,L) = Sn 2 k 2 £ £ K<& n ( K ) j 1 - [(l - |)kp]} dKdz ' (2- 11 ) 

where p is the distance from the optical axis, k = In/A is the optical wavenumber, J 0 
is a Bessel function of the first kind, order zero, and ®„(/c) is the atmospheric index 
of refraction spectrum model being used. For analytic tractability, the Kolmogorov 
atmospheric spectrum is assumed for the following definitions. 

By substituting the Kolmogorov spectrum into Eq. (2.11), further substituting that into 
Eq. (2.10), and solving for the e~ x point one can find the spatial coherence radius of the 
field, p 0 . However, the more commonly used term, r 0 = 2.1p 0 , called the seeing parameter, 
can be shown to be [5] 


r 0 ,sph 


0.423^(zf 



( 2 . 12 ) 
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This is also referred to as Fried’s parameter after being introduced by D. L. Fried [9]. 
With respect to AO, r 0 is a critical parameter for designing a system. Apertures larger in 
diameter than r 0 require AO to see additional gains in imaging resolution. Further, the 
spatial density of wavefront estimation and correction by an AO system are often related to 
r 0 . To be effective, roughly one WFS subaperture and DM actuator spacing should span r 0 . 

To describe the temporal coherence of the field aberrations, Greenwood analyzed an 
ideal AO system with a single pole control filter. He calculated the characteristic frequency 
of the system with respect to the ratio of aperture diameter to r 0 [19] and derived what is 
now called the Greenwood frequency: 

r r L i 3/5 

f g = 0.254fc 6/5 J C 2 n {z)v 5l \z)dz\ , (2.13) 

where v(z) is the perpendicular wind velocity along the propagation path. 

2.1.2 Wavefront Sensing. 

The first step to mitigate the atmospheric effects on light is to be able to estimate the 
distortions introduced for a given propagation path and instant in time. This process is 
called wavefront sensing and requires some form of reference beacon transmitting light 
through the turbulence column of interest. With longer electromagnetic wavelengths, such 
as in radar, the fields of interest can be directly measured and compared. Unfortunately, in 
optics the short wavelengths and corresponding high frequencies are beyond the detection 
range of current electronics. In optics, wavefront estimation requires inferring information 
about the received beacon field from other physical processes [41]. The following 
subsections give an overview of three such methods. 

2.1.2.1 Tilt Tracker. 

One of the most simplistic wavefront sensors is a tilt sensor. Shown in Fig. 2.3, this 
sensor is meant to estimate the low-order wavefront aberrations of x and y tilt. By focusing 
incoming light onto a quadcell detector or sensor array, the deflection of the focused light 
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can be calculated by the centroid equations 




. E ;=1 yiy _ 

Ay= i~^- y - 


(2.14) 


where x and y denote detector column and row positions, x and y are the mean column 
and row, and I x and I y are the summed detected irradiance values for each column and row. 
Using the results from Eq. (2.14), the incoming wavefront’s mean tilt can be determined as 


e x = Ax/f, d y = Ay/f, (2.15) 

where / is the focal length of the focusing optic used, shown in Fig. 2.3. 

Tilt comprises approximately 87% of all atmospheric aberrations [27], and as such, 
using a dedicated fast-steering mirror (FSM) to compensate tilt alone offloads a large stroke 
requirement from the high-order wavefront corrector [25]. Further, the tracker can be used 
with broadband light, has no moving parts, and is computationally simple to operate. These 
reasons make the tracker a very robust wavefront sensor; however, it can only sense the 
lowest order wavefront aberrations. 


2.1.2.2 Shack-Hartmann Wavefront Sensor. 

The Shack-Hartmann WFS uses a tiled array of trackers across the aperture. This 
divides the full aperture into subapertures, whose incident wavefronts are assumed to 
be roughly planar. With this architecture, local wavefront tilts, or gradients, are sensed 
which are then pieced together using a wavefront reconstruction algorithm - discussed in 
Sec. 2.1.3. Because of this aperture averaging effect with respect to high order wavefront 
aberrations, it is advisable to have at least one subaperture per r 0 spacing of the wavefront 
being sensed [25]. 

Figure 2.4 shows a Shack-Hartmann design using square lenslets that focus down on to 
2x2 pixel regions on a larger focal plane array. This is a typical design, but other variations 
exist for a range of applications. One common design decision is to include a band of 
unsensed pixels between subapertures to prevent light from one subaperture ‘spilling’ into 
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Figure 2.3: Principles of operation for a tilt sensor where 9 is the incoming wavefront mean 
tilt with respect to the sensor’s optical axis, d is the lens diameter, / is the effective focal 
length, and Ax is the lateral shift of the focal pattern irradiance. 


an adjacent subaperture’s sensed pixels. However, this requires the use of a larger and more 
expensive focal plane array and typically increases the required sensor readout time. This 
wavefront sensor is very widely used as it shares the same benefits as the tracker does. 
The primary downside of this WFS is its high sensitivity to alignment errors. As in most 
WFS’s there is a tradeoff to be made between dynamic range of detectable wavefront tilts, 
measurement accuracy, signal sensitivity, and speed. 

2.1.2.3 Self-Referencing Interferometer Wavefront Sensor. 

The Self-Referencing Interferometer (SRI) WFS uses interference effects to directly 
sense an incoming field, as opposed to most WFS designs which measure some form of 
wavefront gradients. It is based on the point diffraction interferometer (PDI), shown in 
Fig. 2.5. The PDI focuses all incoming light to a partially transmissive filter in the focal 
plane. The on-axis point in the filter is completely transmissive, acting as a spatial filter. 
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Figure 2.4: Diagram of the Shack-Hartmann WFS. 


When the light is recollimated, an interferogram between the aberrated input wavefront and 
spatially filtered reference wavefront can be seen [24,34], 

Unfortunately, the single interferogram produced by a PDI does not give a direct 
means of forming a wavefront estimate, but this is where the SRI extends. Here, the 
incoming light is split into two separate paths before being recombined. The first path 
is responsible for spatially filtering the incoming wavefronts, similar to the pinhole filter 
used in the PDI. However, in the SRI a single-mode fiber is used instead of a pinhole filter. 
Differently from the PDI, this filtered reference beam is further split up into separate beams, 
and each beam receives a different phase shift - explained below. The second path applies 
no filtering or altering of the beam, but is carefully matched in length to the reference 
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Figure 2.5: Diagram of the Point Diffraction Interferometer. 


legs to maintain a high degree of temporal coherence. Finally, each phase shifted beam is 
recombined with the beacon wavefront to form an interferogram for each. The irradiance 
pattern of each interferogram is measured on an imaging sensor for use in calculating the 
wavefront phase. Figure 2.6 shows a simple layout for an SRI. 

A direct measurement of the input wavefront phase can be calculated using the 
separate interferograms [11,34], Assuming a perfectly coherent source, the m-th sensed 
pixel of the n-th interferogram follows 



In(x m ,ym) =11 \ u b(x,y) + U r (x,y)e Wn \ 2 dxdy. 


(2.16) 


where Ub is the input beacon field, U r is the reference field, and A m is the area of the 
corresponding pixel. If we assume the amplitudes of both fields are equal to A/2, Eq. (2.16) 
can be shown to equal 


In(x m ,y m ) = A [l + COS {(f>bx m ,y m + #n)] ■ 


(2.17) 


Here, the planar nature of the reference wavefront is clear, while the aberrated input 
wavefront remains a function over x m and y m . There are multiple methods for calculating 
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Figure 2.6: The Self-Referencing Interferometer. 


the input wavefront phase using different numbers of interferograms with different 
reference phase shifts applied [11]. For this research, a four bin approach is used. In 
this configuration, phase shifts of 0, n/2, jt, and 3n/2 are applied respectively to each 
reference beam, each interferogram constituting a bin. Following Eq. (2.17), the measured 
irradiances of each interferogram become 


h(x m ,y m ) 

= 

A+A cos (<p b ( Xm ,y m )) 

h(x m ,y m ) 

= A[ 1 + COS (<f> b ( Xm ,y m ) + |) 

\= A-A sin ((f>b( Xm ,y m )) 

h(x m ,y m ) 

= A [l + COS (<f> b (x m ,y m) + 7r)] 

= A- A cos (<p h ( Xm , ym )) 

h(x m ,y m ) 

= A [^1 + COS (<p b (x m ,y m ) + y) 

] = A + A sin ((f>b( Xm ,y m )) 


From this system of equations, the aberrated input wavefront can be easily calculated 
according to Eq. (2.19) and (2.20): 


tan(0i (Wm) ) = 


sin(^ (w „,)) 

COS ((f> b ( Xm ,y m) ) 


WxmJm) ~ h{x m ,y m ) 
h(x m ,y m ) ~ h(x m ,y m ) 


(2.19) 
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( 2 . 20 ) 


fan = Tan -1 • 

For some applications, including the strong turbulence reconstruction techniques in 
Sec. 2.3, it is important to note that the complete field estimate is available with the 
relationship 

U S R, = (h-h) + i(h ~h), (2.21) 

from which the field irradiance can be calculated. 

Equation (2.20) reveals a remarkable point of interest for real-time AO: using an 
SRI, the wavefront phase can be directly measured by merely subtracting the measured 
interferograms, one division, and performing an inverse tangent operation on that result. 
Further, for a spatial shifting SRI these operations can be performed independently for 
each subaperture, meaning that calculations can begin immediately as pixel data from each 
interferogram is read in from the sensor(s). To capitalize on this advantage, thought must 
be given to the pixel readout order of the sensor(s). The computational time required by 
the arctangent operation is a tiny fraction of the full wavefront reconstruction required for a 
Shack-Hartmann. However, the output of the arctangent is wrapped [-/r, n\, which is well 
suited for a segmented DM, but is difficult to reproduce over a continuous face-sheet DM. In 
the case of a continuous face-sheet DM, one needs to implement a wavefront reconstruction 
algorithm in the control system [10,14,30]. Due to the influence functions of a continuous 
face DM, sharp wrapping cuts are difficult to reproduce and may even cause damage to 
hardware. Unfortunately, the reconstruction process eliminates the efficiency gains seen 
over the Shack-Hartmann WFS. 

The SRI still has one major advantage over other WFS technologies in that it is nearly 
immune to the effects of scintillation seen in strong atmospheric turbulence. A scintillated 
field has null points in irradiance over the phase wavefronts that are being estimated by a 
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wavefront sensor. At these points of zero irradiance field estimation issues arise in other 
wavefront sensing technologies [14]. For the SRI, there may not be any beacon light over 
a given subaperture, but the reference beam will have a fraction of the incoming aperture 
averaged intensity. Therefore the local phase is directly calculated as 

* - Tan " (^)= Tan “ = Tan “ G)= 0 (2 - 22) 

In this sense, the SRI is seen to be nearly immune to scintillation; however issues still 
arise in coupling light into the reference fiber when these conditions exist. To minimize 
coupling inefficiencies, an effective tilt tracking and correction system is essential when 
using an SRI. Previous research also investigated using doped reference fibers to achieve 
an optical gain for the reference beam path. Previous research by Rhoadarmer found that 
just enough gain to overcome sensor readout and quantization noise is ideal [35], but in 
practice this is not often implemented. This research seeks to exploit the benefits of the SRI 
with respect to scintillation in conjunction with a novel wavefront reconstruction technique, 
discussed later in Sec. 2.3. 

2.1.3 Wavefront Reconstruction. 

It is often necessary to further process the WFS outputs to achieve an unwrapped 
wavefront estimate. This process is called wavefront reconstruction. For the SRI, WFS 
phase outputs are first sheared into slope space. Other WFS designs typically directly 
output wavefront gradients. As slopes, each measurement point is related to the others 
through a set of linear equations - called the geometry matrix. This system of equations is 
then solved to find the continuous phase wavefront that minimizes the squared error of the 
original system of equations. The following subsections describe this process in detail. 


19 



2.1.3.1 Linear Algebra. 

The problem of relating wavefront gradient measurements to unwrapped phases can 
be written as, 

Ac = b, (2.23) 

where b is the measured gradients in column vector form of length ra, x is an n-length 
column vector of the unwrapped phases being estimated, and A is a geometry matrix 
relating the two of size rax n. The specific form of A will be discussed in Sec. 2.1.3.2, 
but it is a known system parameter. The process of solving for x varies depending on the 
specific properties of A. If m = n and A is full rank, then the solution is straight forward as: 


A- 1 Ax = A~ l b 
x = A~ l b. 


(2.24) 


However, in AO this is rarely the case. Instead, there are traditionally more measurements 
than phase points being estimated (ra > n). This set of linear equations is overdetermined, 
requiring an approach to reduce the parameter space to a best fit solution. Premultiplying by 
A T , (A T A) becomes invertible, giving a best fit solution in terms of the minimized squared 
error for jc [40]: 

A T Ax = A T b 


(A T A)~\A T A)x = (A T A)~ 1 A T b 


x = (A T A)~ l A T b 


x = A'b 


(2.25) 


This is the basic process which is used in the following subsections to estimate wavefront 
phase from gradient measurements from the WFS. Other methods exist for solving a system 
of linear equations such as this; however, the Moore-Penrose pseudoinverse (often referred 
to as the least squares pseudoinverse) is the most commonly used algorithm. 
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2.1.3.2 Zonal Reconstruction. 


Zonal reconstruction uses a simple relative relationship between measured slopes and 
the phase estimation points. This relationship is captured by the Geometry matrix, of which 
there are two predominant types depending on the way in which a particular WFS senses 
wavefront gradients. The Fried geometry is typically associated with the Shack-Hartmann 
WFS, while the Hudgin geometry is used with the lateral shearing interferometer (LSI) and 
SRI. 

Wavefront sensors implementing the Fried geometry sense wavefront slopes at the 
center of each subaperture and relate those slopes to phase estimation points at the corners 
of each subaperture through Eq. (2.26) and (2.27) [12]. Figure 2.7 shows a 2x2 subaperture 
array in the Fried geometry. 

( S X -Fried)m,n = ~^{4 > tn,n+\ + 0m+l,n+l) — (fpmji + 4 > m+\,n) (2.26) 

(S y -Fried)m,n = T^m+hn + </>m+l,n+l) ~ ^(<P»un + </>m,n+ 1) (2.27) 

In this geometry, there are two slope measurements per subaperture, creating 2 Q equations 
for P phase estimation points. Relating this to the discussion in Sec. 2.1.3.1, this system of 
linear equations can be written as 


G f </> = 




(2.28) 


Alternatively, wavefront sensors implementing the Hudgin geometry estimate phases at 
the center of subapertures based on gradient measurements between adjacent subaper¬ 
tures [21]. Figure 2.8 demonstrates a 2x2 subaperture wavefront system. The mathematical 
relationship for this geometry is more simplistic than the Fried geometry as, 


($x-Hudgin)m,n ~ fim.n+l tym.m 

(2.29) 

( Sy-Hudgin)m,n = <Pm +1 ,n ~ tym.n- 

(2.30) 
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Similarly to the Fried geometry, this creates 2 Q equations for P phase estimations, creating 
the system of linear equations: 


G H (f> = 


S x-Hudgin 
Sy-Hudgin 

Both Eq. (2.28) and (2.31) can be solved as in Sec. 2.1.3.1 by: 


0 = G^s, 


(2.31) 


(2.32) 


where G and s are the appropriate geometry matrix and slopes from above. 

Further, G f is usually precomputed for the given wavefront sensing geometry being 
implemented. This eliminates the need to compute the pseudoinverse at each frame. The 
resulting matrix is called the reconstruction matrix, 

R = G\ (2.33) 


This poses the general least squares wavefront reconstruction problem: 

4> = Rs. (2.34) 


2.1.4 Wavefront Correction. 

The process of correcting the aberrated field inputs is heavily driven by the application 
of interest. Techniques exist for correction of both field amplitude and phase [2,7,36]; 
however, only phase compensation is discussed further here. Traditional AO uses the 
methods described in the previous sections to estimate the phase wavefronts and eliminate 
any estimated aberrations. Figure 2.9 depicts the physical process of phase conjugation. 
With an accurate estimation of the aberrated input, a conjugate shape is applied to the 
surface of a wavefront corrector such that when the light reflects, the wavefront is near 
flat. As mentioned previously, the low order aberrations of tilt are often corrected 
separately from the higher order aberrations using a combination of deformable mirror 
and FSM [25,41], 
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Figure 2.7: The Fried Geometry. Slopes are measured at the center of the subaperture and 
are related to the phases at each comer. Slope vectors shown here include both x and y 
components for each subaperture. 


There are a number of kinds of wavefront correction devices that lend to specific 
applications [41]. This research uses a two inch diameter Optics In Motion FSM [3] 
and 12 x 12 actuator Boston Micromachines Micro-Electro-Mechanical System (MEMS) 
DM [1], 

2.1.5 AO Control. 

To implement a practical AO system, some form of control scheme is necessary to 
derive useful commands to the wavefront correction device(s) based on the estimated input 
aberrations [25,41]. In practice, most AO systems rely on closed-loop control for the 
increased robustness when integrating complex beam control systems. Figure 2.10 shows 
a high-level control diagram for a tracking system. An aberrated beam entering the AO 
system will have some time-varying wavefront tilt, <p t . The beam first reflects off of the 
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Figure 2.8: The Hudgin Geometry. Slopes are measured between adjacent subapertures, 
which are related to the phases at the center of each subaperture. 


FSM, which leaves the beam with some residual tilt error, (p t -res (near zero for a good 
control system in steady-state operation). The residual tilt error is sensed and estimated by 
the tracker. The estimated error is then fed into a control law to calculate the next set of 
adjustments to the FSM - with the ultimate goal of driving the sensed error to zero. For the 
AO system used in this research, a first-order linear filter of the form, 

c k = ae k -\ + fick-u (2.35) 

is implemented in both the tracking and higher-order control systems. In Eq. (2.35), c k is 
the controller output at time step k, e k is the sensed error at time step k, a adjusts the DC- 
gain, and /? is a servo leakage factor. In a pure integrator setup, /? = 1, but slightly lowering 
this value decreases sensitivity to any unsensed modes or misregistration errors. 

Following the control law, a digital-to-analog converter (DAC) takes the digitized 
controller calculations and converts them to a usable analog voltage. This is typically 
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Figure 2.9: Wavefront phase conjugation using a deformable mirror. 


buffered and/or amplified prior to being sent to the wavefront corrector. The FSM being 
used in this research implements its own internal closed-loop controller to maintain the 
mirror position to the analog position commands [3]. Based on the position of the FSM 
with respect to the optical axis, a tilt is imparted to the incoming beam, and the process 
repeats for each frame in time. 

2.1.6 Metrics of Performance. 

The holy grail of almost any optical system design is to achieve diffraction limited 
imaging performance for the application at hand. For this to be the case, there must be zero 
wavefront aberrations throughout the optical system - including the imaging path. With 
this in mind, the goal of AO is to remove any observable wavefront aberrations from the 
light of interest before reaching the science sensors or target. In AO, the primary measure 
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Figure 2.10: High-level control scheme for a tilt compensation system. 


of correction performance is Strehl ratio. Strehl ratio is defined as [38] 




(2.36) 


where ( H(f x ,f y ) is the aberrated optical transfer function (OTF) [18] and H d [(f x ,f y ) is the 
diffraction limited OTF. Both of these OTF’s are in terms of their spatial frequencies f x and 
f y . There are several approximations to this form that better suit various applications or 
measurement techniques [25,38]. For this research field-estimated Strehl is used. 

As the name suggests, field estimated Strehl takes a given field and computes an 
approximate Strehl without the need to propagate to a focal plane. This is often useful 
for efficient numerical simulations or direct Strehl estimation from WFS reconstructions. 
It is given by 

3= I , (2.37) 

Aff \U(x,y)\ 2 dxdy 

where U(x,y ) is the field within an aperture of area A. For the simulations to follow, 
Eq. (2.37) is discretized to the form 


(2 .38) 

mean[\U(x,y)\ 2 \ 

For hardware performance evaluation, field estimated Strehl is subject to the accuracy of the 
wavefront reconstruction. Therefore, when evaluating the effect of different reconstructors 
on AO performance in hardware, field estimated Strehl is not a good metric. An alternative 
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is to use either the tracking camera measurements or a separate, higher resolution, science 
camera to compute r H(f x ,f y ) from the sensed irradiance point spread functions. 

2.2 Scintillation and Adaptive Optics 

In electromagnetics, when a field interferes with itself, areas of both constructive and 
destructive interference result. This is considered a scintillated field [48]. Figures 2.11 and 2.12 
show the irradiance of a point source beacon propagated through both weak turbulence and 
strong turbulence, respectively. In each case, D/r 0 « 2 and A = 1,55/nn, but there are key 
differences between the two that should be highlighted: 


Weak Turbulence Example 

• Cl = 5 x 10“ 15 , L = 1km 

• 0-2 = 0.01 

• Smooth phase function 

• Relatively evenly spread irradiance 


Strong Turbulence Example 

• Cl = 2.2 x 10“ 17 , L = 244km 

• 0-2 = 1.06 

• Phase function contains sharp cuts 

• Scintillated irradiance with both peaks 
and areas of near zero power. 


When describing a given optical path through the atmosphere, weak atmospheric 
turbulence generally implies that any transmitted light will only accumulate a randomized 
phase as a result of any refractive or diffractive effects. Figure 2.13 shows light being 
spherically emitted from a star that is roughly planar upon reaching Earth’s atmosphere. 
However, once these planar wavefronts propagate through the turbulent atmosphere, the 
wavefronts become aberrated. Strong atmospheric turbulence generally refers to an optical 
path that causes not only phase accumulation but also interference effects. For this to 
happen, the beam accumulates a phase aberration early in the optical path. With further 
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(a) (b) 


Figure 2.11: Recieved point source irradiance(a) and phase(b) propagated through weak 
turbulence. 



Figure 2.12: Recieved point source irradiance(a) and phase(b) propagated through strong 
turbulence. 
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Figure 2.13: Effects of the randomized atmosphere on planar phase wavefronts. 


propagation the aberrated wavefronts begin to constructively and destructively interfere - 
causing scintillation in the received light. 

Unfortunately for AO, scintillation presents a unique set of problems with respect to 
phase conjugation. Fried wrote the first significant publication about scintillation’s effect on 
the corresponding phase function [13], and he found that in areas of no light within the field 
a new phase component forms. Such points are called branch points. Figure 2.14 shows 
an idealized field that has a single branch point at the center. At the branch point there is 
zero amplitude, and the phase goes undefined. To make the field continuous in a modulo- 
one wave sense, a rotational phase emanates outward from the branch point, producing 
the corkscrew shape. Expanding to the general case, Fried showed that all branch points 
create a rotational phase field which is in addition to the generally associated irrotational 
wavefront aberrations from propagations through the turbulent atmosphere [13]. This can 
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be written simply as 


t = (f> irr + trot, (2.39) 

where trot is the rotational component of the phase due to branch point effects and tirr 
is the irrotational component typically associated with atmospheric optical propagations. 
Note that in weak turbulence, <f>rot = 0 because there are generally no branch point effects 
present. 

As seen in the theoretical branch point of Fig. 2.14, the rotational phase is defined 
over a one wave region. This, combined with the rotational nature, requires a sharp cut 
emanating outward from the branch point along the wrapping edge of the one wave window. 
This is called a branch cut and can be seen in parts of the reconstructed phase in Fig. 2.12. 
Branch cuts always link pairs of branch points or travel from one branch point to outside 
the pupil region. The rotational phases are defined on an arbitrary one wave region, and the 
2n window can be shifted without changing the modulo-27r congruence of the wavefront. 
However, doing so moves any branch cuts. For example, the rotational phase of Fig. 2.14 
is windowed between -n to n. If the window is changed to 0 to 2n, the wrapping boundary 
would move to the opposite side of the corkscrew, and the branch cut would shift sides 
of the field. Note that shifting this window by integer multiples of 2n give equivalent 
placements of the branch cuts. Likewise, if the wrapping window is kept fixed, a uniform 
piston added to the entire field has the same effect [31,45]. This latter technique is preferred 
for this research simply due to computational efficiency and can be written as 

t' rot = W [trot + h ], (2.40) 

where t' rot is still modulo-2^- equivalent to trot■ Figure 2.15 gives a graphical example of 
this process. 

Unfortunately for most AO, the wavefront corrector used is a continuous face-sheet 
design. These correctors will inevitably have a fitting error along any branch cuts in the 
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Figure 2.14: Theoretical rotational phase function surrounding a centered branch point in 
a propagated field. 


field being conjugated. The following section poses a strategy to minimize the impact of 
branch cut fitting errors from the wavefront corrector. 

2.3 Wavefront Reconstruction in Strong Turbulence 

As described in Sec. 2.1.3, raw WFS measurements usually need reconstruction prior 
to correction. Herein lies a major problem for traditional AO techniques when branch 
points are present in the observed field. Any reconstruction matrix calculated using the 
least squares pseudoinverse from Eq. (2.25) cannot reconstruct a rotational component [13]. 
This is not a function of the particular WFS/DM geometry chosen - it is instead a behavior 
of the Moore-Penrose pseudoinverse. Figure 2.16 demonstrates this behavior. 

Fortunately, this is not a dead end. There are actually a number of algorithms one can 
use to recover (f > rot . Pellizzari provides the most complete comparison to date with a focus 
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Figure 2.15: Effects of adding a uniform piston, h, to the rotational field while maintaining 
a fixed wrapping window. Modulo-2^- congruence is maintained, but the branch cut moves 
as h changes. 


on the AO application [31]. Of the many algorithms from a range of science applications, 
the primary constraint that must be considered for use in AO is computational efficiency. 
The wavefront must be reconstructed in real-time at the frame rate of the WFS - often in 
the kilo-frames per second. Pellizzari used this to narrow down his evaluation to fifteen 
reconstructor options for which he did a more thorough comparison. The research that 
follows builds from his conclusions. In the subsections below, the most promising real¬ 
time branch point-tolerant wavefront reconstructors are described. 

2.3.1 Post-Processing Congruence Operation Reconstruction. 

The most promising class of algorithms is based on a post-processing congruence 
operation (PCO) to find and add back in the rotational portion of the field [31,45], This 
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Figure 2.16: Applying a least 
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reconstructor to a rotational field. 


research chose the SRI because it directly senses the wavefront phase - including any 
rotational components present. We can thus formulate the problem in terms of Eq. (2.39) 
and the PCO technique as 

0rrr = 0LS = RlS T ( </>SRI ) > 

and (2.41) 

4>rot = 'W[<f>SRI-<l>irr\, 

where <Psri is the SRI’s raw phase measurement and the 170) operator finds the x and 
y gradients of 0, consistent with the reconstructor geometry being used [13, 15,31,45]. 
Substituting Eq. (2.41) into Eq. (2.39) gives 

0PCO = 0LS + [05- <I>ls] ■ (2.42) 

This process maintains the traditional least squares phase reconstruction technique for the 
irrotational phase component and finds the modulo-2^- equivalent rotational phase out of 
the original measurement. This process adds only a very small computational overhead to 
gain the rotational component reconstruction. 

There is further room for improvement on the basic PCO process. As discussed in 
Sec. 2.2, the branch cut locations can be varied by applying a piston to the rotational field 
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prior to wrapping Eq. (2.40). The subsections that follow take the general PCO algorithm 
and attempt to optimize the particular realization of rotational phase with respect to the 
branch cut locations by seeking an optimal piston addition, h. 

To compare realizations of rotational phase, a cost function posed by Yenema [45] 
and refined by Pellizzari [31] called intensity weighted cut length (IWCL) is used. 
IWCL first locates any branch cuts in a rotational field-again by shearing the raw phase 
to get wavefront gradients, and searching for sharp gradients. It then integrates the 
SRI’s estimated field irradiance from Eq. (2.21) along the estimated cuts. Results from 
both Yenema and Pellizzari show high negative correlation between IWCL and Strehl 
ratio [31,45]. This research deviates slightly from the previous IWCL implementations 
in not normalizing the integrated IWCL by the total field irradiance. While the normalized 
IWCL is useful for comparing specific test cases, there is no advantage to the additional 
computations in a real-time AO system. 

By seeking to minimize IWCL, branch cuts are moved towards areas of minimized 
irradiance, thus keeping the effects of DM fitting errors along the cuts to a minimum. 
This is only one cost metric - others can be chosen based on the application of interest. 
One alternate was developed for beam projection applications, which is presented later in 
Sec. 4.6. 

23.1.1 LSPV. 

Least squares principle value (LSPV) was first published with respect to AO by 
Venema and Schmidt [46]. LSPV+n takes the one wave range of possible h offsets, and 
divides it n times: 

>W = [o:^:2jr). (2.43) 

While the definition above splits h evai equally across the one wave range, there are variants 
that do not necessarily follow this constraint. The LSPV+n reconstruction process then 
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evaluates all n values of h eva i in terms of IWCL to determine the best (f> mt : 

<f>LSPV+n = (f>LS + [ (f>SRI ~ <pLS + h opt ^ , (2.44) 

hop, = argmin /w [IWCL (<f> rot (h eval ))\ . (2.45) 

LSPV+1 and LSPV+4 were the primary algorithms investigated in this level of PCO 
optimization. LSPV+1 may seem like a poor choice, but with only one value of h eva i, there 
is no need to compute IWCL - giving a computational advantage. LSPV+200, or other 
large n, can be used to test for the true ideal h opt for a given field but is not suggested for 
real-time implementation due to the computational requirements and diminishing marginal 
gains. Figure 2.17 shows one particular LSPV+4 execution. It is obvious that some values 
of hevai are less optimal than others. 

2.3.1.2 Iterative Probing. 

Proposed by Pellizzari [31], this algorithm starts a narrow LSPV+5 search about the 
midpoint between the previous frame’s h opUk - \ and (h opt ). In this case, the narrow LSPV+5 
search evaluates 

7 [ 1 h op tjc~ i + ( h op t ) 

hevaii = [ -0.3 -0.15 0 0.15 0.3 J +-^- waves. (2.46) 

The calculated h optk from the above LSPV+5 search serves as a seed for a second search. 
The new search narrows the range of h eva i for a finer investigation of the IWCL parameter 
space about h opUk \ 

hevaii = [ -0.07 -0.04 0 0.04 0.07 ] + K pUk waves. (2.47) 

Pellizzari motivated this by showing that the IWCL-space is roughly parabolic. This 
iterative probing method allows for further precision with still minimal additional 
computational requirements. 
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(C) (d) 


Figure 2.17: These plots depict the rotational portion of a wavefront in radians as a piston 
shift is applied [45], Branch points are marked by X’s and O’s based on their respective 
charge, and branch cuts are highlighted with white. The non-unique nature of the rotational 
phase is highlighted by applying piston: n (a), nil (b), 0 (c), and -n/1 (d). Plots adapted 
from [45]. 


2.3.1.3 Histogram Search. 

Also proposed by Pellizzari [31], this algorithm minimizes the area of (f> ro , near the 
wrapping boundaries. It first calculates a histogram of the </> rot phase values. From there, it 
calculates h' opt as the shift requried to place the minimally sized histogram bin at the nearest 
wrapping edge. Additionally, a narrow LSPV+3 search explores the IWCL parameter- 
space about h' opt to arrive at the final h opt . 
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2.3.1.4 Brent’s Method. 


The final PCO algorithm posed by Pellizzari utilizes what is known as a Brent’s 
Method [33]. Unlike the iterative probing algorithms above, this type of search algorithm 
minimizes the required evaluations of the cost function for the optimization problem at 
hand. For this, it means fewer evaluations of IWCL through a golden ratio search of the 
h parameter-space. Knowing that the IWCL parameter-space is roughly parabolic, this 
algorithm recursively refines a parabolic fit by adjusting h. Each parabolic fit uses four 
data points, from which the search algorithm divides that search range in half. Successive 
searches maintain a ratio between cost function evaluation inputs such that previous IWCL 
evaluations are reused as the search narrows. 
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III. High-Fidelity Simulation Development 


This chapter highlights the development of a Matlab model of AFIT’s own AO testbed. 
From the start of development, this was intended to be a high-realism wave-optics based 
simulation of the true system dynamics. The following sections highlight the development 
decisions, lessons learned, and configuration capabilities for the simulation model. 

The simulations were developed in the Matlab programming environment. As this is 
a discrete numerical model, careful consideration was given to adequate spatial sampling 
of the electromagnetic fields involved. Although this discussion does not cover the theory 
behind numerical simulation of wave-optics, this information can be found in Schmidt [38]. 
Additionally, Matlab toolboxes, WaveProp and AOTools, were used to augment the 
methods proposed by Schmidt. These toolboxes were developed by the Optical Sciences 
Company (tOSC). Much self-authored functionality and integration code was developed 
for this research and is described in the sections that follow. 

3.1 AFIT Adaptive Optics System 

AFIT’s AO system testbed was built as a multi-purpose platform on which to give a 
hands-on experience for Air Force students being introduced to AO. As it has evolved, it has 
become an extremely capable and configurable platform for advanced AO research [37], 
Despite this, having a simulation model on which to conduct initial research, where many 
variables of experimental work are idealized, creates a powerful time and money saving 
tool. The following subsections introduce the physical system and how it was translated 
into the computer model. 

3.1.1 Hardware Specifications and Setup. 

Built by SAIC, AFIT’s AO system consists of two independent control systems. As 
described in Sec. 2.1.2, a tracking system handles the low order tilt wavefront aberrations, 
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while a separate high-order system senses and corrects the remaining wavefront errors. 
Figure 3.1 shows the optical layout of the AO system which can be fed by anything that 
relays to an entrance pupil at the FSM. As is currently setup, the beam diameter at the FSM 
is 14mm before being reduced by a factor of three in a 4-/ relay to the DM. From then 
on, the beam remains at 4.67mm and is relayed once more to the sensor entrance pupil 
which precedes the Shack-Hartman WFS beam-splitter. This size gives ten DM actuators 
across the beam as shown in Fig. 3.2. A border of actuators remain that are slaved to 
the actively driven ‘master’ actuators. Slaving unsensed actuators minimizes influence of 
adjacent actuators at the beam edges on the actively driven actuators [43], For this research, 
the Shack-Hartmann WFS is being completely bypassed and going unused. Instead, the 
higher-order WFS is a 19 x 19 subaperture spatial SRI operating at 1.9Kfps with one frame 
of latency from sensor exposure to DM commands out. The tracking sensor has a field 
of view of 1.8mrad over a 32 x 32 sensed pixel array operating at 4.76Kfps, again with 
roughly one frame of delay from exposure to FSM command out. Note that this field of 
view is referenced to the AO system entrance pupil, not the demagnified beam entering the 
actual tracker. 

This spatially-oversampled WFS-DM geometry is important in properly reconstruct¬ 
ing the wavefront with a phase-only WFS when there is potential for large phase differences 
between WFS subapertures [6]. Such large phase differences are a natural consequence 
of branch cuts within the field. Therefore, a high-resolution wavefront reconstruction is 
formed from the 19 x 19 SRI measurements, which is then downsampled to the DM res¬ 
olution prior to the control-law calculations. It is important when downsampling only to 
include phase estimates within the active aperture, otherwise the outer edge of downsam¬ 
pled points will always be artificially biased towards zero. 

Simulating the hardware system allows for a few simplifications. First of all, there is 
no need to simulate any optical relays of the entrance pupil. Instead, all phase corrections 
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Figure 3.2: Geometry of SRI subapertures (dotted squares) to master DM actuators (blue 
dots) the beam extent is shown in red. 
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are performed on the same field, which is then appropriately scaled in irradiance before 
being sensed by the tracker and SRI. Of course, this assumes a perfectly aligned and zero 
aberration optical system; however, if one wanted to add either effect, it would be a simple 
addition. 

In simulation, a one meter telescope collects light from the optical propagations. This 
is not as critical as D/r 0 is, where D is the entrance aperture diameter. As discussed before, 
the ratio of wavefront corrector size to Fried’s parameter, r 0 , is an important consideration 
when designing a system. For this simulation work, any atmospheric conditions can be 
generated, leaving the choice of entrance aperture size somewhat arbitrary. The entrance 
telescope is simulated to have a focal length at the simulated target distance for each 
scenario, giving a collimated output beam to the AO system. 

The telescope also demagnifies the one meter entrance pupil to a 4.67mm exit pupil. 
This design decision ignores the one-third demagnification from the FSM to DM in the 
physical system but increases the simulation simplicity. Due to the changed angular 
magnification at the simulated FSM, the only change between the simulations and physical 
system is the final gain coefficient between the tracker controller to the analog FSM 
voltages. Since no inertial effects of the FSM are simulated, this gain discrepancy has no 
effect on the fidelity of results. To justify excluding simulating inertial effects of the FSM, 
a comparison between the FSM command frequency spectra (at the highest Greenwood 
frequency tested) and the physical FSM’s measured frequency response showed only a 
small content approaching the FSM’s roll-off frequency. 

With the optical system design handled, consideration of the separate control systems 
must be paid. Often in computer simulations of AO systems, the tracking system and high- 
order system operate at integer frequency multiples of one another (i.e. Tracking at 4Kfps, 
High-order at 2Kfps); however, this is not often the case for physical systems. This is an 
issue when calculating the optical fields that feed the system. As these are computationally 
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very expensive to generate, a minimum number of frames should be sought. The solution 
used in this research was to calculate separate lists of sampling times for each control 
system based on their respective sampling rates. With those, a unique merge of the two 
lists is performed, while maintaining a table that links each point in time to which sensor(s) 
is active. This resulting list is not uniformly spaced in time, but it does not have to be. 
Knowing the turbulence wind speeds and point in time, WaveProp is able to generate the 
propagated fields without issues. The listing below gives pseudo code through this process 
for additional clarification. 
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1 % Load propagated fields 

2 % Setup sensors and correctors 

4 for currFrame = ItnumFrames 

5 % Apply simulated telescope 

6 % Apply last FSM and DM corrections 

8 if SRI active 

9 % Detect, reconstruct, downsample, control law, filter 

10 end 

11 if TRK active 

12 % Detect, threshold, centroid, control law 

13 end 

u % Calculate Strehl 

i 6 end 

With the above system coded up, initial simulations began. Unfortunately, the first set 
of Strehl plots looked like the left subplot of Fig. 3.3. As the figure shows, initially the AO 
system is completely off, and the Strehl is roughly zero. The tracking system turns on 0.01s, 
and the Strehl picks up to roughly 0.1. At 0.02s, the high-order system turns on, leading 
to a good initial improvement, but the system quickly goes unstable. Upon investigation of 
the DM commands, it was obvious that tilt was dominating the DM corrections. 

While not apparently an issue in the previous work by Pellizzari where the separate 
control systems operated at the same frequencies, here the higher-order and tilt control 
systems create a condition where they fight each other. As there is always going to be a 
slight residual correction error from the tracking system, the SRI detects this and tries to 
compensate, although at a slower sampling rate. The faster tracking system will have time 
to sense and correct for the same residual that the higher-order system is still processing. 
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When the control loop gains are high enough, this behavior can quickly overwhelm the 
limited stroke of the DM and even drive the whole AO system unstable. 

The solution is simple: ensure there are no common correction modes between any 
control loops. In this case, project out tilt from the higher-order system. In traditional AO, it 
is co mm on to augment the traditional least-squares reconstructor matrix with a tilt-removal 
operation, such that the reconstructed wavefront has no tilt components. The process for 
this is shown in the code below: 


i %% Remove Z—Tilt: 

3 % Generate arbitrary orthogonal tilted wavefronts: 

4 [xt,yt] = meshgrid(l:nDMacts); 

5 % Zero—mean: 

6 xt = (xt - mean (mean (xt) )) ; 

7 yt = (yt - mean (mean (yt) )) ; 

8 % Only consider the active aperture: 

9 xt = xt(activeMask == 1) ; 
o yt = yt(activeMask == 1); 

1 % Transform from 2D-space into Reconstructor-space 

2 % and normalize the tilted wavefronts: 

3 XZ = xt*xt'/(xt'*xt); 

4 YZ = yt *yt'/(yt'*yt); 

5 % Create the Tilt-Removal Operator Matrix: 

6 P = eye(sum(activeMask(:))) - XZ - YZ; 

8 % Combine to form a tilt—removed Reconstructor: 

9 R_tr = P * R; 


44 




Unfortunately, with the PCO-based reconstruction techniques, this process is not as 
simple. As outlined in the previous chapter, the PCO process is 

<j>pco = <Pls + 'dV [<I>sri ~ • (3.1) 

However, to make </>pco tilt-removed, two forms of least-squares reconstruction are needed: 

4>PCO = <f>LS-tr + '3V [(f>SRI - </>LS-ti\ , (3.2) 

where </>Ls-tr is a tilt-removed least-squares reconstruction while <f>Ls-ti retains tilt. This 
does indeed give a tilt-removed reconstruction, but it requires the calculation of two least- 
squares reconstructions rather than only one. 

A better solution exists. Realizing that the control law follows the reconstruction, 
a design decision exists on where to calculate actuator slaving. Often, this is rolled in 
with the least-squares reconstructor; however, this technique will not work with the PCO 
calculation of 4Kot- Additionally, the slaving operation should ideally be the last calculation 
prior to sending the DM commands out. So a tradeoff can be made in this case: only 
calculate a single, tilt-included, least-squares reconstruction which gives a tilt-included 
PCO reconstructed wavefront. Downsample the reconstruction, send it through the control 
law, and apply one final matrix operator which removes piston, removes tilt, and applies 
slaving. This is the optimal data-flow process for a PCO reconstructor; shown in Fig. 3.4. In 
total, this requires just two matrix multiplies - one for the least-squares reconstruction and 
one for the final filtering. With the tracking and higher-order control systems decoupled, 
performance improved greatly. This is seen in the right plot of Fig. 3.3. 

3.1.2 System Radiometry. 

With the simulation model framework mostly in place, focus could be turned to 
modeling radiometry throughout the optical system. This is crucial if noise effects are 
going to be accurately simulated, as photon shot noise is signal-dependent [17] and other 
noise sources are only meaningful in relation to the signal level. 


45 



Decoupled Tilt 



Figure 3.3: The effects of neglecting tilt decoupling on Strehl ratio. Right plot shows tilt 
properly decoupled over the same input wavefronts for comparison. 
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Figure 3.4: Optimal process flow for a PCO based reconstruction algorithm. 
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3.1.2.1 Radiometric Assumptions. 

In the physical AO system, power can be varied to each signal path using polarized 
beam splitters preceded by adjustable half-wave-plates [37]. This requires the input light 
be initially filtered by a polarization analyzer. For the simulation work, no polarization 
effects are modeled - only the power splitting ratios. In this sense, the input light for the 
simulations would be equivalent to just after the initial polarization analyzer on the physical 
table. 

A number of assumptions were made in this radiometric model: 

1. No losses due to reflections. 

2. No losses due to dispersion or misalignments. 

3. Ideal beam splitters. 

Realistically there is a small loss through the system despite the anti-reflective coatings 
used on all optics. Additionally, the beam splitters have a finite extinction ratio. This is 
relevant for the beam paths that theoretically receive no light. Realistically, the Shack- 
Hartmann and auxiliary beam paths receive a small percentage of the light. 

Additionally, realistic sensor models were created in WaveProp based on the camera 
datasheets from Goodrich. For both the tracker and SRI, a Goodrich SU320-KTSW sensor 
is used, with the key specifications shown in Table 3.1. 

3.1.2.2 Realistic SRI Reference Fiber Coupling. 

Previous simulations assumed the SRI had a constant fiber coupling coefficient into 
its single-mode reference fiber, regardless of the input wavefront condition. This is a 
rough approximation for quick evaluation but is insufficient for a realistic noise analysis, 
determining low signal performance, and examining transient behavior before the AO 
system reaches steady-state. 
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Table 3.1: Goodrich SU320-KTSW Sensor Specifications 


Resolution 

320 x 256 f 

Pixel Pitch 

25jum 

Quantum Efficiency 

>0.65^ 

Full Well e“ 

750,000e- 

RMS Read Noise 

183e -ttt 


f Can address a smaller rectangular region of the array for increased speed, 
ft Actual quantum efficiency at 1550nm is likely greater, but 0.65 was used for these 
simulations. A sensor calibration can experimentally determine this value [22], 
ttt Based on correspondence with Goodrich technicians for a FPA gain of one. 

Based on previous work by Wheeler [47], frame-by-frame fiber coupling efficiency 
and exit mode shape can be calculated. Using the physical optic and fiber parameters of 
AFIT’s SRI hardware, the single mode accepted into the fiber (LP 0 i) is back-propagated to 
the SRI entrance pupil. With this, the input aberrated beam is projected onto the known 
back-propagated LP 0 i fiber mode, giving a realistic frame-by-frame coupling coefficient. 
Similarly, the exiting LP 0 i fiber mode is a magnified version of the input mode, based on 
the specific re-collimating optics used after the output of the fiber. This output mode is 
scaled by the current frame’s coupling coefficient as 

U fiberout (k) = (u input (k), U 0l - in ) U 01 - out , (3.3) 

where (•, •) denotes an inner product, k represents the current time step, U input is the input 
aberrated beam of the SRI reference leg at time step k, f/ 0 i-,„ is the normalized back- 
propagated LPqi mode that is coupled into the fiber, and U 0 is a normalized version 
of the collimated fiber exit mode. This operation captures the stochastic variations in 
coupling efficiency, while other effects like Fresnel reflections are treated separately and 
deterministically. As is shown in Figs. 3.5 and 3.6, fiber coupling is particularly important 
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in modeling transient behavior immediately after closing the loop in strong turbulent, low 
signal situations. 

3.1.2.3 Power Splitting. 

As shot noise is dependent on the signal level at the detector and the other noise 
sources are only meaningful relative to the signal level, careful work had to be taken to 
ensure proper radiometry throughout the simulated AO system. This problem was initially 
approached by adjusting the power splitting between tracker and SRI to achieve roughly 
equivalent signal levels on each detector. For the tracker this is straight forward, but the 
SRI presents a challenge with the input-dependent fiber coupling. For this calibration, 
a low-Rytov turbulence scenario was used to determine the optimum tradeoff in power 
splitting. Figures 3.5 and 3.6 illustrate this tradeoff. Before the AO system is engaged, 
the fiber coupling is near zero. First the tilt system is engaged, giving a very slight 
increase in coupled light. This coupled signal must be sufficient to form the SRI wavefront 
estimate, otherwise the higher-order control loop may not successfully close when engaged. 
When the higher-order system engages, coupling immediately jumps as the DM starts 
correcting. This is shown in Fig. 3.5 when the fiber power becomes stronger than the 
signal power. However, as fiber coupling efficiency is directly correlated with Strehl ratio, 
stronger turbulence degrades the power coupled into the SRI reference fiber, as seen in 
the medium-Rytov case of Fig. 3.6. A balance must be struck between supplying enough 
power to initially engage the higher-order control loop while not taking too much power 
from the signal leg. The power splitting can thus be optimized for the expected atmospheric 
conditions, or as in this research, a general compromise was found that accommodates a 
wide range of scenarios. Figure 3.7 shows the power splitting breakdown through the 
AO system that was empirically determined, and as it turns out matches up closely to the 
physical AFIT hardware setup. 
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Time [s] 


Figure 3.5: SRI reference fiber coupling as the AO system starts up. The left plot shows 
the power in the fields for both the signal and reference (labeled Fiber Power) legs through 
the SRI. The right plot calculates the fiber coupling efficiency over time. Atmospheric 
conditions are f g = 20Hz, Rytov number = 0.04. 






Figure 3.6: SRI power splitting between reference and signal legs on the left using a 
LSPV+4 reconstructor for turbulence with a Rytov number of 0.488 and Greenwood 
frequency of 20Hz. The resulting fiber coupling is shown on the right. 
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Collected Light 



Figure 3.7: Approximate power splitting within the simulated AO system. The SRI’s 
interferograms are represented with /, through / 4 . I r and I s are the reference and signal 
paths, respectively. The fiber loss will depend on the instantaneous Strehl ratio of the field 
entering the SRI. 
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3.1.3 Sensor Noise Effects. 

With the system radiometry setup, focus can be turned to the noise sources. For this 
discussion, focus is on a single pixel of one SRI interferogram, as we are most interested 
in noise effects relating to the reconstruction performance. However, it should be noted 
that the balanced power distributions and same detector models used in the SRI and tracker 
yield roughly equivalent SNR’s. Noise sources analyzed include photon shot noise <r] hot , 
detector readout noise cr 2 read , and ADC quantization noise cr 2 quantization so that the total system 
noise is given by 

= y/^hot + a lad + Quantization' ( 3 ‘ 4 ) 


Photon shot noise follows a Poisson distribution and is only a function of the signal 
strength. Its noise contribution to Eq. (3.4) can be expressed as [35] 

Ihot = K 


ffqe_ 

N hin: 



(Is + Vch) 
2 


(3.5) 


where K is the mean photoelectron count for a given pixel, rj qe is the detector quantum 
efficiency, N hins is the number of SRI interferograms captured simultaneously, T int is the 
detector integration time in seconds, Jl px is the area of one pixel in square meters, r/ c is the 
mean coupling efficiency into the single-mode fiber, and I s and I r are the irradiances for the 
signal and reference legs of the SRI, respectively. As this is an interferogram, I s and /,- are 
averaged, rather than simply adding their signals. 

As mentioned in Table 3.1, 


0- 2 read = ( 183 e f . 


(3.6) 


This read noise number varies depending on the gain and calibrations used, but this gives 
a rough middle ground for an InGaAs camera in this class and represents one bit of digital 
resolution. This sensor uses a 12bit ADC to quantize the analog photocurrents. Assuming 
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a rounding ADC, the quantization noise can be expressed as [35] 


Quantization = T^LSB 2 

12 (3.7) 

= ^(183e-) 2 , 

where there are 183 photoelectrons per least significant bit (LSB) of the 12bit ADC. 

With all of the noise sources defined, the challenge becomes to define a meaningful 
SNR to test AO and reconstructor performance. The noise is now defined in terms of 
photoelectrons for an individual pixel of one interferogram. The signal will thus be defined 
as the average number of photoelectrons incident on one pixel over the same integration 
period. Making use of Fig. 3.7, roughly 5% of the entrance power reaches the pixel from 
both the signal and reference legs of the SRI. Since this is an interference pattern, the 
combined signal could range between 0-10% with a mean about 5%. Therefore, the signal 
estimate is 5% of the input power. Note that this decreases if Strehl performance degrades 
for any reason. Putting Eqs. (3.4)-(3.7) together with the previous definition of signal, SNR 
can be defined as 

SNR = = (3.8) 

y^shot ^read quantization 

SNR = K — . (3.9) 

Using the above equation for SNR, the required input beam power can be calculated 
for 1,55/rm light. Note that these SNR values are not in decibels. In the extreme, before 
the higher-order loop is closed, actual signal levels can be much lower than these estimates 
due to poor fiber coupling, as seen in Figs. 3.5 and 3.6. Because of this, initially closing the 
loop has been found to be the primary obstacle for the SRI-based AO system to overcome. 
Ensuring the tracking system is closed before trying to close the higher-order loop is 
essential in low signal conditions. 
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3.2 Atmospheric Propagations 


While not the primary focus of this research, realistic aberrated optical inputs are 
needed to test the AO system over. An ideal point source beacon at a wavelength of 1,55/rm 
was used in all cases. As mentioned before, care was taken in ensuring adequate spatial 
sampling of the wavefronts using methods from Schmidt [38]. WaveProp also alerts users if 
inadequate sampling is detected. Table 3.2 lists the spatial sampling used in this research, 
where cr 2 is the Rytov number, d/r 0 is the ratio of the DM actuator spacing to Fried’s 
parameter, 8 src is the beacon source plane grid spacing, 8AOsys is the receive plane grid 
spacing entering the simulated telescope, N Grid Points is the number of discrete spatial 
samples per side of the simulated field, C 2 is the index of refraction structure constant along 
the propagation path, and D z is the propagation path length. 

3.2.1 Parameter Exploration. 

To fully characterize the AO system performance, a wide parameter-space was 
investigated. Most importantly, this research evaluates the effect of an increasing Rytov 
number on reconstruction performance. Additionally, for each Rytov number, a range of 
Greenwood frequencies and SNR’s are investigated. Finally, two r 0 ’s were considered at 
d/r 0 = 1 (Challenging) and d/r 0 =1/2 (More realistic). This is all summarized in Fig. 3.8. 


Table 3.2: Atmospheric propagation setup parameters. 



d/r 0 

8 src [mm] 

$AOsys [mm] 

N Grid Points 

C 2 n 

D z [km] 

0.04 

1 

0.4 

8 

1024 

2e-14 

1 

0.5 

1 

4 

4 

2048 

le-15 

20 

1 

1 

5 

7 

2048 

4e-16 

50 

0.5 

0.5 

9 

13 

1024 

7.88e-17 

80 

1 

0.5 

12 

15 

2048 

2.2e-17 

244 
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• Horizontal Propagation of a point source 

• Hold d /r 0 constant in all cases: 1, ¥z 
■ Vary Rytov Number: 0.04, 0.5, 1 

• Vary Greenwood Frequency: 20,40, 60, 80, 100Hz 

Lead to various combinations of D z , C£, and V w i nd 


Kvind 



Figure 3.8: Point source atmospheric propagation setup. 
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IV. Simulation Results 


This chapter first presents a compiled set of simulation results from the various 
reconstruction methods introduced in Ch. II. Obvious performance improvements from the 
traditional least squares reconstruction approach are seen, especially as turbulence strength 
increases. While that was the original intended scope for the simulation portion of this 
research, further investigations led the author to several additional insights on PCO-based 
reconstruction that greatly clarify previous questions and hypotheses. These new findings 
include least squares phase anomalies due to gradient aliasing, the cause and potential 
solutions to sharp fades when using PCO reconstruction in strong turbulence, and refined 
behavior models of h opt during closed-loop operation. This chapter concludes with an 
extension of the PCO technique in strong turbulence to beam projection applications. 

4.1 Reconstructor Performance 

Due to the large number of simulations run, this section presents a summarized set 
of the raw results. While exploring the effects of SNR and Greenwood frequency were 
informative along the research pathway, the more general performance trends between 
each reconstruction technique can be explored without these details. The tables that 
follow present averaged statistics of the field-estimated Strehl ratio from each one second 
simulation run. For each Rytov number and reconstructor combination, all simulations 
were grouped together to calculate the mean, standard deviation, and normalized standard 
deviation of the field-estimated Strehl ratio. The mean Strehl ratio measures absolute 
AO correction performance, and the standard deviation of the Strehl ratio measures the 
stability of the AO correction. Higher mean Strehl ratios and lower standard deviations are 
desired. More detailed plots of the simulations run, linking effects of SNR and Greenwood 
frequency, are found in Appendix A. 
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Table 4.1: Averaged field-estimated Strehl ratios over all Greenwood frequencies and 
SNR’s tested for the cases listed. 


<4 

LS 

LSPV+1 

LSPV+4 

Histogram 

Iter. Probe 

Brent’s 

0.04 

0.5406 

0.5416 

0.5416 

0.5416 

0.5416 

0.5416 

0.5 

0.3242 

0.3844 

0.3869 

0.3876 

0.3872 

0.3846 

1 

0.0949 

0.2047 

0.2177 

0.2183 

0.2202 

0.2134 


Table 4.2: Field-estimated Strehl ratio standard deviation, averaged over all Greenwood 
frequencies and SNR’s tested for the cases listed. 



LS 

LSPV+1 

LSPV+4 

Histogram 

Iter. Probe 

Brent’s 

0.04 

0.0329 

0.0322 

0.0322 

0.0322 

0.0322 

0.0322 

0.5 

0.0835 

0.0498 

0.0496 

0.0470 

0.0487 

0.0488 

1 

0.0677 

0.0647 

0.0674 

0.0625 

0.0626 

0.0646 


Table 4.3: Field-estimated Strehl ratio normalized standard deviation, averaged over all 
Greenwood frequencies and SNR’s tested for the cases listed. 



LS 

LSPV+1 

LSPV+4 

Histogram 

Iter. Probe 

Brent’s 

0.04 

0.0628 

0.0613 

0.0613 

0.0613 

0.0613 

0.0613 

0.5 

0.2581 

0.1310 

0.1295 

0.1223 

0.1267 

0.1278 

1 

0.7162 

0.3198 

0.3152 

0.2899 

0.2881 

0.3058 
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4.2 Low SNR Reconstruction and Phase Discrepancies 

Among the first observations made was in the initial weak turbulence testing. In 
weak turbulence, there are no branch point effects present, and therefore <p rot - 0 in the 
absence of any noise effects. The expectation was to observe equivalent performance 
between all reconstructors. However, a consistent performance improvement of all tested 
PCO reconstructors was seen over the least squares reconstructor, specifically at the lowest 
SNR’s tested. This was very much a surprise, as prior work hypothesized that the addition 
of noise may be debilitating to the PCO technique [31]. This sparked a series of follow-up 
investigations to explore what was going on. All tests performed verified the initial findings 
that all the PCO reconstructors outperformed the least squares reconstructors in low-SNR 
situations. The question then became: Why? 

An initial theory was proposed by the author that was based on differing noise 
propagation through the two reconstruction techniques. The hypothesis is based on the 
idea that the least squares operator could actually amplify the SRI measurement noise such 
that 

\\<f>L Sn \\ 2 > m\ 2 (4.i) 

where </> LSn is the least squares propagated noise component, (f) n is the noise component 
from the measurement process, and || • || 2 is the Frobenius norm [16]. This idea was formed 
by calculating the condition number [16] of the least squares reconstruction matrix and 
finding it to be roughly 15 - ignoring the suppressed reconstruction modes. Now relate this 
to the PCO algorithm defined in Eq. (2.42): 

0PCO = *pLS + *41? [(/>SRI ~ </>ls ] • (4.2) 

The above form can be expanded in terms of all noise components involved: 

frpCO = + $LS„ + 4V {TV [<f> + <f> n ] - (^LS^ + 0LS„)} • (4.3) 
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The theory proposes that the two <f> LSn terms cancel each other out in the PCO process, 
leaving only <p n contributing to the final reconstruction. Under the assumption of Eq. (4.1), 
only having <p n would be less of a degradation than A rough analytical trace for 

this theory was created, but it needed additional verification to justify some simplifications 
made. 

A simulation was written to model this process which proved the theory wrong. This 
verification simulation models AO system behavior in steady-state closed-loop operation 
by generating an irrotational test field (could be as simple as a plane wave) that has 
a nearly flat wavefront. The mean field irradiance is scaled using the equations from 
Sec. 3.1.3 to have an SNR of roughly 1.5. This field is sent through the SRI model to 
add realistic measurement noise effects, and that field, </>sri, is sent through both a least 
squares reconstructor and a simple LSPV+1 reconstructor. From these, the noise terms can 
be extracted and examined. 

Specifically, the simulation found that 

<f>LS n = <f>n, (4.4) 

which invalidated the author’s original theory for the PCO reconstructor performance 
improvement. In hopes of gleaning more, additional realizations of noise were run, 
and something unexpected was noted in a number of cases. Amongst other things, the 
simulation plotted <f> LSn - <p n , usually showing no difference. However, Figs. 4.1a and 4.1b 
show two examples of what appeared to be branch points at first glance. 

Yet, these are not branch points for a number of reasons. First, the field that was input 
to the SRI had a uniform irradiance. One could argue that shot noise creates areas of zero 
irradiance, which is a valid argument. The detected field irradiance can be calculated from 
the SRI interferograms by Eq. (2.21), and Fig. 4.2 shows the estimated input irradiance 
corresponding to the measurement of Fig. 4. lb. It is clear that areas of zero detected 
photons are scattered all over the aperture, but there is no correlation to the anomalies 
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observed in Fig. 4.1b. Since these effects do not seem to be caused by the input irradiance, 
looking to the SRI phase measurement, a clue is found. Shown in Figs. 4.3a and 4.3b, the 
points from Figs. 4.1a and 4.1b can be traced to points in equivalent locations within the 
SRI phase that make transitions larger than n radians between adjacent subapertures due 
to measurement noise. While examining <(>sri, closer examination was given to both the 
least squares and LSPV+1 reconstructed wavefronts, and it was clear that while the PCO 
reconstructor successfully reconstructed the original phase, the least squares reconstruction 
had anomalies at the locations of the sharp phase differences in </>sri- 

To further examine this phenomenon, a second simulation was run that isolates the 
noise generated large phase differences from the SRI. A hypothetical <f> S Ri is generated that 
is completely flat except for two adjacent center subapertures that are set to -0.9 n and 0.9n, 
as shown in Fig. 4.4a. Without the surrounding additional noise, this simulation exposes the 
aperture-wide reconstruction effects due to large noise-induced adjacent phase difference. 
Figure 4.4c shows the LSPV+1 reconstruction, matching the SRI wavefront with an extra 
n piston shift from the PCO h term. On the other hand, the least squares reconstruction 
clearly identifies the source of the phase anomalies evident in Fig. 4.1. 

The (f> ro , term from the PCO algorithm cancels out the least squares anomaly; however, 
this exposes a misnomer in this case. The wrapped term in the PCO algorithm, which in a 
noiseless case is equivalent to any rotational component in the phase function, can actually 
be any phase components missed, or in this case artificially introduced by the least squares 
reconstruction relative the SRI phase measurements. The author likens this to accounting 
for slope discrepancies when using a Shack-Hartmann WFS. Since the SRI measures phase, 
not wavefront gradients, the author terms this as accounting for the phase discrepancy. In 
the case of the noise induced phase discrepancies, they are not actually rotational in nature. 
As described in Sec. 2.3 a least squares generated operator cannot reconstruct rotational 
phases. Further, Fig. 4.5 plots the gradients for one of these anomalies, where red represents 
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a positive field and blue is negative. The author stresses these arguments to say these are 
not so-called “noise-induced branch points”. They are not rotational, and they are actually 
anomalies introduced by the least squares reconstruction that do not exist in the true phase 
function. 

Since these phase discrepancies result from noise creating artificially large adjacent 
phase measurements on an otherwise nearly-flat wavefront, it is now clear why all PCO 
algorithms outperform the least squares reconstructor when no branch points are present 
and only at the very lowest SNR’s. The author’s original theory failed to explain why the 
performance gain was only seen at SNR’s lower than five, as it predicted a lowered noise 
propagations for any noise present. The new theory requires such low SNR’s to induce 
the large adjacent phase measurement errors. If only using a least squares reconstructor in 
these cases, the control system will attempt to conjugate these phase discrepancies that do 
not exist in the true wavefront leading to potentially unstable AO performance. 

The cause for the observed least squares behavior turns out to be an aliasing effect 
of the wavefront gradients fed to the the least squares reconstructor. Section 4.5 presents 
the research that led to this recent understanding. Any gradients larger than n radians are 
assumed to be wrapping cuts; however these points are not wrapping cuts. So, the large 
adjacent gradients are wrapped, but the surrounding gradients remain. This is the cause for 
the phase discrepancy surrounding such points. 


4.3 Closed-Loop Behavior of h opt 

Pellizzari’s work found a numerical PDF approximation for h opt [31], which found 
iKpt) = 0. However, there is a very important difference between preceding research and 
this work in the implementation of the PCO wrapping operator. Previously, the wrapping 
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<pLS„ ~ <t>n 


<t>LS n ~ 0r> 



5 10 15 20 25 30 35 5 10 15 20 25 30 35 


(a) (b) 

Figure 4.1: Different realizations of unexpected differences between LS noise term and 
PCO noise term. 



corresponds to the measurement of Fig. 4. lb. 


62 



















<PSRI 

m r 

5 % 

10 . 1 

1 

1 

10 

<PSRI 

15 . - 

■ 

| 

15 

1 0 tg, 


20 . 


1 £ 20 


35 "J b ™k] ■ 

1 

_ _1 25 

1-2 30 

la 5 

1 - ■ . /-ijj 
' ■ 


5 10 15 20 25 30 35 5 10 15 20 25 30 35 


(a) (b) 

Figure 4.3: SRI phase measurements of a tilted wavefront where SNR« 1.5. 
subplots, (a) corresponds to Fig. 4.1a and (b) corresponds to Fig. 4.1b. 
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I: 


In the 


was performed by 


^wm = 


Tan 1 {imag [exp (i2n(f >)\, real [exp (i2n(f>)]} 
In 


(4.5) 


which is physically correct but not computationally efficient. Switching to a modulus 
operator, this process can be implemented easily on most CPU’s and FPGA’s used for 
AO control. This is simply written as 


TV [0] = mod ((f). In). 


(4.6) 


Both methods are correct, but the wrapping window for Eq. (4.6) is now from 0-2n radians, 
rather than -7T-7T radians. When this research began testing LSPV+1, based on Eq. (4.6), the 
AO performance was terrible. After recalculating the PDF of h opt , it was obvious that the 
PDF calculated by Pellizzari needs to be defined with respect to the wrapping window being 
used. By shifting his entire PDF to be over the new wrapping window, (h opt ) = n. It turns 
out the initial implementation of FSPV+1 tested here, used the least optimal h, leading 
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Hypothetical <p SRI 





Figure 4.4: Simulation of low-SNR noise effects between least squares and LSPV+1 
reconstructors. 



Figure 4.5: Simplified wavefront gradient diagram of a noised induced phase anomaly. The 
center two squares represent the central discrete pixels of the anomaly. 
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to the poor results observed. This is an important lesson for practically implementing 
these algorithms or anyone implementing a non-standard wrapping function for a PCO 
application. 

Over the course of this research, it was noted that the LSPV+1 reconstructor 
demonstrated near equivalent AO performance to the more advanced PCO reconstructors 
tested. The computational simplicity of LSPV+1 makes it a promising lead for high-speed 
and low-latency AO. It was also noted that the LSPV +1 performance gap between the other 
PCO algorithms widens as the Rytov number increases. Although the performance is still 
very similar, as seen by Table 4.3, the increasing performance gap between LSPV+1 and 
the other PCO reconstructors at a 2 x = 1 prompted an investigation as to why. 

LSPV+1 always uses h = (h opt ), which in the weak Rytov case is more than 
sufficient when no branch points are present. However, as the Rytov number increases and 
additional rotational phase components are introduced, one can surmise that only having 
one realization of <f> rot could be limiting in some cases. The previous research did not 
show significant variations in the PDF of h opt as the Rytov number increased above weak 
turbulence, perhaps due to the absence of noise effects. However, a series of simulations 
were run to retest this idea with the hypothesis that the PDF of h opt broadened away from 
(h opt ) as the Rytov number increased. To generate the PDF’s, LSPV+200 was run over 
a series of one-second cases, each with increasing challenge for the AO system. For 
each simulation, a histogram was generated of the h opt values chosen by the LSPV+200 
algorithm to generate a rough PDF. Figures 4.6 and 4.7 show the simulation results with 
the histograms as well as field estimated Strehl ratio plots overlaid with the normalized 
h opt as a function of time. As the AO beacon conditions worsen, the histogram shows a 
broadening of hop ,’s PDF. All turbulence parameters showed a slight broadening effect of 
the distribution, but increasing Rytov number is by far the dominant factor. 
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The confirmed broadening of the h opt PDF explains why the LSPV+l-LSPV+4 
performance gap widens as the atmospheric conditions worsen. This fact also may lend 
to choosing a specific PCO reconstructor based on expected conditions. As the conditions 
become increasingly difficult, particularly with respect to Rytov number, it may be more 
beneficial to take the additional computational hit with one of the other PCO reconstructors. 

4.4 Fades in Strong Turbulence 

While debugging the reconstruction and AO system simulation, a mechanism was 
implemented to record videos of many important system parameters as the system stepped 
through time. Following the development phase, this capability was used to study the 
overall behavior of the AO system when using a PCO wavefront reconstructor. This 
additional work resulted in identifying the source of fades when using a PCO reconstructor 
when a significant rotational phase is present in the input wavefronts. 

At system startup, the WFS observes the fully aberrated beacon wavefronts, <f>. At 
this point, the complete rotational phase from any branch points reaches the WFS. When 
the AO system switches on, the DM conjugates the reconstructed WFS field. Now, the 
WFS senses only the residual wavefront error, <f> res id - not the fully aberrated wavefronts. 
The control law maintains the DM commands at roughly the conjugate of the incoming 
aberrations, but the wavefront reconstructor is no longer seeing the fully aberrated field. 
As <l> r esid flattens out, the PDF of h opt becomes increasingly centered about (h opt ). 

This was all known prior to the research presented here; however, additional insights 
may be gleaned by thinking deeper into the process. In steady-state, closed-loop operation, 
the residual wavefront error is near-zero. Any branch points and corresponding branch 
cuts that are present within 0 may not be visible in <f> resi d, and even if they are partially 
apparent, they would most likely be of reduced visibility due to partial compensation. In 
this regime, the branch cut finding algorithm may easily miss cuts that are actually present 
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crjj =0.488, d/ro =0.5, SNR=5 



Figure 4.6: LSPV+200 Strehl ratio and h opt plots over one second simulation runs to 
evaluate the behavior of h opt '% PDF as atmospheric conditions change. 
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g\ =0.488, d/r 0 =1, SNR=60 



Time [s] h Offset [waves] 







Figure 4.7: LSPV+200 Strehl ratio and h opt plots over one second simulation runs to 
evaluate the behavior of h opt '% PDF as atmospheric conditions change. 
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in the corrected field. The author poses a question as to what IWCL really means, since 
IWCL was meant to be calculated over the fully aberrated wavefront, <p, not <f> res id- 

So what is going on when a PCO algorithm changes h op , while in closed-loop 
operation? Well, it depends primarily on the rotational effects present in the field; in other 
words, the Rytov number. Starting the discussion with a weak Rytov number case, the 
rotational component of the PCO reconstruction was seen to be almost always completely 
zero, with the exception of any noise induced effects presented in Sec. 4.2. In this weak 
turbulence case, h opt very rarely changes, and when it does change due to noise induced 
phase anomalies, it has no effect on the overall system performance. 

In the medium Rytov number cases simulated, similar behavior is noted. Figure 4.8 
shows a plot of field-estimated Strehl ratio, overlaid with each frame’s h opt for an LSPV+4 
reconstructor. In this regime, h opt changes frequently, but overall AO performance remains 
unaffected. This initially diagreed with the author’s theories, as the PCO rotational 
phase component was seen to have consistent notable content. However, the most recent 
discoveries show that this content is not rotational in nature. The cause is discussed at 
length in Sec. 4.5. Without a significant rotational component within the residual phase 
error, h opt varies in a way to keep the PCO phase discrepancies within the wrapping window 
and away from the wrapping boundaries. At a change in h opt , only the reconstructed 
wavefront’s piston changes, but since piston is removed prior to sending DM commands 
out, there is no performance impact. While there are a few exceptions, this is the primary 
behavior observed in the medium Rytov number cases. 

Things change with strong Rytov number cases. As Fig. 4.9 shows, any changes in 
h opt mostly result in sharp fades. More specifically, the fades distinctly follow changes in 
h opt . Using the video capability developed for debugging, it was noted that the tracking spot 
shape is not effected at these fades but is shifted off-axis. This behavior was consistently 
observed at the fades in Strehl ratio. It seemed that the change in h opt somehow propagates 
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Time [s] 


Figure 4.8: The top plot shows field estimated Strehl ratio and h opt over time with an 
LSPV+4 reconstructor. The bottom plot shows the IWCL for each of the four <p rot 
realizations over the same time interval. This is in conditions of a 2 = 0.5, f g =40Hz, 
d/r 0 = 1/2. 
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Time [s] 

Figure 4.9: The top plot shows field estimated Strehl ratio and h opt over time with an 
LSPV+4 reconstructor. The bottom plot shows the IWCL for each of the four <p rot 
realizations over the same time interval. This is in conditions of cr^ = 1 ,f g =40Hz, 
d/r 0 = 1/2. 


a change in tilt through the AO system, thus causing the fades. This prompted several 
simulation tests, one of which lead to an explanation. 

A number of AO simulations were rerun with the LSPV+4 reconstructor, this time 
calculating the tilt components for each PCO rotational phase realization. Since h opt 
seemed to be at the root of the behavior, the least squares reconstruction was ignored for 
these tilt calculations. Figure 4.10 shows the results from the three Rytov cases. First to 
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note is that the tilt components of the PCO rotational phase definitely differ as h is varied! 
For the weak case, all rotational tilt components vanish as soon as the higher-order control 
loop closes - as expected. In the medium case, the three most likely h-v alues produce 
small, yet similar tilts. However, in the strong case, even after the higher-order control loop 
engages, there is widely varying tilt components between each rotational phase realization. 
This is the key. 

In the strong Rytov case, when the PCO algorithm changes h opt , a new wavefront 
tilt realization is reconstructed. When this propagates through the rest of the higher-order 
control data path, described in Sec. 3.1.1, the final tilt-removal operator cuts out whatever 
tilt component exists. To the tracking system, these changing tilts are certainly measurable. 
This is the cause for the fades in Strehl ratio. While the higher-order correction of the 
wavefronts do not suffer during the changes in h opt , it takes the tracking system roughly two 
frames to account for the tilt changes from the high-order system. To explicitly demonstrate 
this process, Appendix B traces a hypothetical rotational phase through the AO system. 

4.4.1 Potential Fixes. 

All the tilt components are known variables, and as such, a solution to prevent the 
fades is certainly possible. In fact, there are likely a number of potential solutions. One 
method would be to implement a feed-forward tilt-offloading path to the tracker from the 
higher-order system. Another method is to allow the higher-order control system to correct 
for some or all of the reconstructed tilt. At first, the author recalled the tilt coupling issues 
presented in Sec. 3.1.1, but by lowering the control system gains, a stable state can be 
reached. While this imposes larger stroke requirements on the deformable mirror, it is a 
relatively simple controls problem compared to the feed-forward implementation with two 
different operating frequencies. 

The AFIT AO system model was modified to compare AO performance over 
equivalent realizations of high Rytov turbulence with both a tilt coupled and decoupled 
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Figure 4.10: Tilt magnitudes for each PCO rotational phase realization of a LSPV+4 
reconstructor over time for weak, medium, and strong Rytov numbers. The higher-order 
control system is seen to engage at t =0.008s. 
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control system. The limitations of such a setup are made clear by selecting three 
increasingly aggressive control gains. Figure 4.11 shows the results, where the left column 
is decoupled and the right column is coupled. The gain starts at relatively conservative 
points, giving good performance in each case. However, note there are no fades in the 
coupled case, even as h opt changes! Increasing the gain, the decoupled performance 
sharpens but the coupled performance begins to go unstable. Further increasing the gain, 
the decoupled system still runs well, but the coupled system quickly drives itself unstable. 
Note that this test in simulation allows for perfect component alignments. In reality there 
will always be a slight misalignment between components and sensors. Because of this, 
a zero DC gain controller would be required on both the tracker and higher-order control 
system to implement a coupled tilt control setup. 

This test verifies the theory on tilt propagation from the changing PCO rotational phase 
realizations to the tracking system as the cause for the sharp fades previously observed with 
these reconstructors. Further, suggestions for alleviating this are given. 

4.5 High SNR Phase Discrepancies 

A very recent modification to the simulations was to separate and plot <f> LS and <p rot 
as the AO frames increment. To the author’s great surprise, phase discrepancies similar 
to those described in Sec. 4.2 were clearly apparent in high SNR conditions of all Rytov 
numbers tested. This section describes the author’s initial tests and theories as to why 
these phase discrepancies exist and what impact they have on AO performance. The author 
stresses that these are only initial theories, and further work is needed to refine these ideas. 

Figure 4.12 shows an example of the AO system in weak turbulence, with an SNR 
of 60, prior to closing any of the correction loops. In such low Rytov number and high 
SNR conditions the PCO rotational phase term was expected to have no content, but that 
was not the case. The first clue came when the control loops started correcting in the weak 
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Figure 4.11: Plots of AO performance with tilt decoupled (left column) and tilt coupled 
(right column) between the two control systems. Each row implements an increasing 
higher-order control gain. 
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turbulence case. Within several frames, the PCO rotational component became flat - as we 
had expected for the entire low-Rytov simulation, shown in Fig 4.13. 

The primary differences with regard to the reconstruction process between the two 
above cases are the magnitudes of the residual wavefront gradients. These phase effects 
are also present in the higher Rytov number cases, but the low Rytov case excludes any 
rotational field effects from the analysis. A number of simplified tests were conducted 
to look at how high localized wavefront gradients propagate through the reconstruction 
process. The first question was to explore how the SRI and reconstructors handle slopes of 
increasing magnitude. As previously described, the slopes that the reconstructors operate 
on are generated by shearing the raw SRI measurements in the x and y directions. To 
account for wrapping cuts, the sheared slopes are wrapped to be within -n < s < n radians. 
This slope wrapping operation sets up the fundamental limit on slope magnitudes that 
can be properly reconstructed by the least squares process. Any slopes between adjacent 
subapertures larger than n radians will be interpreted, incorrectly, as a wrapping cut. 

Figure 4.14 explores this idea with two tilted wavefronts in the first row of plots. 
The corresponding SRI measurements of each tilted field are shown in the middle row, 
and the final row shows each least squares reconstruction. The first column is correctly 
reconstructed, but the second column, which exceeds the tilt sampling limit, has two 
significant errors. First, the magnitude of tilt over the wavefront is greatly reduced from its 
actual tilt. Second, the direction of tilt is reversed. The author hypothesizes that the degree 
of each of these errors will change based on how much tilt is being passed into the system 
- similar to aliasing in discrete signal theory [29]. 

Figure 4.15 extends the previous aliasing test closer to the atmospheric simulation 
conditions. A flat wavefront with a sharp, linear slope down the center of the field 
is created. The slope is at roughly O.Sn radians per-SRI subaperture. Additionally, 
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Figure 4.12: Plots of AO state in conditions of cr 2 = 0.04, d/r 0 = 1. Both tracking and 
high-order control loops are open. 
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Residual Phase Error SRI Phase Measurement 



Figure 4.13: Plots of AO state in conditions of cr 2 x = 0.04, d/r 0 = 1. Both tracking and 


high-order control loops are closed. 
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Hypothetical Test Phase Hypothetical Test Phase 



Figure 4.14: Inducing aliasing with the SRI and least squares reconstructor. The left 
column shows roughly 0.9 n radians of tilt between adjacent subapertures, while the right 
column has roughly 1 .5n radians of tilt between adjacent subapertures. 
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a central region of the slope is set at roughly 1.5/r radians per subaperture. This is 
similar to most of the atmospheric wavefront being within the slope sampling limit of 
the reconstructor, with small regions that exceed it. The first observation made was that 
the least squares reconstruction mostly resembles the input wavefront, but it is clear there 
are some discrepancies, including aliasing over the high-slope region. In the initial tests, 
only the LSPV+1 algorithm was used for comparison. While it does produce a correct 
modulo-2^- congruent reconstruction, the <f> rot realization created has sharp cuts that are less 
than ideal for any continuous face-sheet DM. This observation motivated a test with the 
LSPV+4 reconstructor, which calculated the most accurate reconstruction. To the PCO 
algorithms, the sharp cuts from the LSPV+1 </> rot realization are treated as if they were 
branch cuts. In the LSPV+4 reconstructor, the cuts are penalized, giving a realization that 
minimizes the sharp transitions. 

The author claims these tests point out that sharp, localized wavefront gradients are 
the cause for the observed phase discrepancies when reconstructing a high-SNR field. The 
LSPV+4 rotational phase term from Fig. 4.15 has a high resemblance to the features seen 
in the atmospheric simulations. While these features were not found in low-Rytov number 
closed loop operation, they were found in the higher Rytov number closed loop simulations. 
Figures 4.16 and 4.17 show closed loop operation for cr 2 = 0.5 and a 2 = 1, respectively. 
The strong turbulence plots show the system as it is initially closing the loop, and some 
rotational phase due to branch points is likely present in addition to the phase anomalies 
being discussed here. 

These tests shed additional light on the increased performance drop between LSPV+1 
and the other PCO reconstructors as the Rytov number reaches one. The broadening 
of h opt ’s PDF must be, in part, due to these phase discrepancies. Further, the previous 
simulations performed by Pellizzari would likely not have encountered these issues with a 
simulated SRI of over twice the subaperture density. While increasing the SRI subaperture 
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Hypothetical Test Phase SRI Phase Measurement 



Figure 4.15: Inducing aliasing with the SRI and least squares reconstructor. A tilt runs 
down the center of the wavefront, and the very central portion exceeds the SRI tilt sampling 
limit. This is traced through the SRI and various reconstructors. 
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density is a potential solution to this behavior, one must consider the additional splitting 
of light and its impact on SNR. The initial results shown here demonstrate an unintended 
benefit of the PCO algorithm as it compensates for reconstruction errors due to localized 
wavefront gradient aliasing. 


4.6 Application to Beam Projection 

All research up to this point regarding PCO-based wavefront reconstructors has been 
studied with respect to one-way AO systems. Such one-way applications include, but are 
not limited to, remote sensing and astronomy. Another class of applications involve dual¬ 
path propagation, namely those which not only receive light, but also project compensated 
outgoing beams. Beam projection applications include optical communication and laser 
weapon systems. These applications typically involve more horizontal propagation paths, 
thus entering the strong turbulence regime much more frequently than vertical path 
applications. This section presents initial results and observations for a range of wavefront 
correction techniques with respect to beam projection applications. 

The author hypothesizes most of the same conclusions made for the single-path AO 
problem hold true but introduces an alternative to intensity weighted cut length (IWCL) 
for PCO reconstruction. The author hypothesizes that cut length (CL) is a better PCO 
reconstruction metric for beam projection applications where the outgoing beam is not 
amplitude modulated to match the received beacon field. While all PCO reconstructed 
fields will be modulo-2^- equivalent, different PCO variations differ in the placement of 
detected branch cuts within the field. Assuming an ideal corrector, this does not actually 
make a difference. However, when a continuous face-sheet DM is used, there will always 
be a DM fitting error in trying to reproduce the sharp branch cuts. 

For initial comparison of wavefront reconstruction techniques with respect to beam 
projection applications, a simplified simulation was developed using tOSC’s WaveProp. 
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Again, an ideal point source beacon illuminates the turbulent path of interest, with a Rytov 
number of one. The receive telescope is 40 centimeters in diameter, focused at the beacon 
distance. An idealistic SRI senses the collimated wavefront (no fiber coupling or noise 
effects), and an ideal DM applies the appropriate conjugation to the outgoing beam. The 
ideal DM has no stroke limit, negating the need for a FSM or tracking system. No temporal 
effects are included in this initial experiment. Instead, unique realizations of turbulent paths 
are generated, through which the beacon propagates. For each turbulent realization, the AO 
system senses the field and various reconstruction techniques calculate corrections for the 
outgoing beam. Each of these pre-compensated beams are then back-propagated through 
the turbulence to the beacon source plane, where the irradiance patterns are averaged over 
all turbulence realizations for each reconstruction method. The author acknowledges the 
lack of temporal effects will not paint a complete picture of expected AO performance. 
However, this is an important first step to both develop a dual-path test simulation, as well 
as compare general reconstructor performance in these difficult atmospheric conditions. 

The following tests compare reconstruction performance of traditional least squares, 
LSPV+20 using IWCL as a metric, LSPV+20 using CL as a metric, and Xphase as 
developed by tOSC. LSPV+20 was chosen instead of LSPV+4 in order to more accurately 
compare the PCO technique with respect to each metric function. Lor reference, ideal phase 
conjugation and ideal field conjugation are also tested. The first round of simulation results 
is shown in Pig. 4.18. To the author’s initial surprise there was no apparent difference 
between the ideal phase conjugation, PCO methods, and Xphase. This points out an 
important point: these reconstruction algorithms all produce modulo-2^- equivalent results, 
but each uses a different rotational realization in hopes of optimizing placement of any 
branch cuts with respect to an expected DM fitting error. This initial simulation did not 
include any DM inter-actuator coupling effects. Purther, the SRI and DM were assumed to 
have the same spatial resolution as the simulated fields themselves. Without any of these 
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realistic effects, all of the above reconstruction methods gave exactly equivalent results. 
The poor performance of the least squares reconstruction and great performance of the 
field conjugations also serve as good sanity checks before moving forward. 

The simulation was modified to include realistic SRI sampling and DM parameters 
(except actuator stroke). As this introduces some subjectivity for the parameters chosen, 
a range was evaluated to explore the current realistic parameter-space. Inter-actuator 
couplings of both 10% and 20% were evaluated for DM’s of size 10 x 10, 19 x 19, and 
32 x 32. The respective SRI subaperture arrays were 19 x 19, 37 x 37, and 63 x 63, 
again, following the sampling recommendation set by Barchers [6]. The WFS to DM 
command downsampling method was a matrix operator based on a pyramid convolution 
three subapertures wide, only considering the active subapertures. The ideal phase and 
field conjugation methods remain at the full field sampling density for comparison. 

Again, 100 realizations of turbulence were run for each of the six DM setups described 
above. The results of which are shown in Figs. 4.19-4.24. Additionally, Figs. 4.25 and 4.26 
provide a visual example of the ‘best’ and ‘worst’ DM fits for an example turbulence 
realization for each reconstructor. These initial results show some interesting behaviors; 
although, not all good. Most notably, the 10 x 10 DM seems to outperform the higher 
density DM’s. Upon reinspection of the simulation code, the AFIT DM setup was used - 
including a guardband ring of DM actuators which is not used in the higher density DM 
setups. While it is not definitive that the guardband was the cause for the apparent improved 
performance, it is a difference that must not be included in future tests. Secondly, the 
LSPV algorithms and Xphase often show higher peak on-axis power than the ideal phase 
conjugation. The author hypothesizes this is a result of a constructive interference effect; 
although, that is only a theory. One final concern upon viewing these results was the varying 
power levels between each simulation case. In the initial code, the turbulence seeds were 
set based on the computer clock - resulting in the observed behavior. For future work, the 
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seeds are now set based on the trial number such that each DM simulation will use the same 


set of turbulence realizations. 

These initial results show the field conjugation performing best and least squares 
worst. Beyond that, Xphase seems to slightly edge out LSPV+20 IWCL in most cases. 
Surprisingly, the LSPV+20 using CL as the reconstruction metric underperformed the 
IWCL metric. The author hypothesizes this is another interference effect at play, but is 
uncertain on the precise means. However, before making too many conclusions based on 
the results at hand, the author stresses the flaws pointed out in the preceding paragraph. 
Further simulations are being conducted to form a better result-set for comparison. 
Additionally, future simulations including temporal effects are crucial to better understand 
this problem. 
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Ideal Field Conjugation 


Ideal Phase Conjugation 





Figure 4.18: Target-plane irradiance patterns, averaged over 100 turbulence realizations 
with an ideal DM. The central circles show the diffraction limited spot size. 




















Ideal Field Conjugation 


Ideal Phase Conjugation 





Figure 4.19: Target-plane irradiance patterns using a 10 x 10 DM with 10% inter-actuator 
coupling. 


















Ideal Field Conjugation 


Ideal Phase Conjugation 
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Figure 4.20: Target-plane irradiance patterns using a 19 x 19 DM with 10% inter-actuator 
coupling. 
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Ideal Phase Conjugation 
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Least Squares Reconstructed 
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Figure 4.21: Target-plane irradiance patterns using a 32 x 32 DM with 10% inter-actuator 
coupling. 
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Figure 4.22: Target-plane irradiance patterns using a 10 x 10 DM with 20% inter-actuator 
coupling. 
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LSPV+20 using IWCL 
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Figure 4.23: Target-plane irradiance patterns using a 19 x 19 DM with 20% inter-actuator 
coupling. 
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Figure 4.24: Target-plane irradiance patterns using a 32 x 32 DM with 20% inter-actuator 
coupling. 
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Figure 4.25: Example DM profiles using a 10 x 10 DM with 20% inter-actuator coupling. 


95 




















-200 -100 0 100 200 (deg) 

(mm) 


-200 -100 0 100 200 (deg) 

(mm) 


Figure 4.26: Example DM profiles using a 32 x 32 DM with 10% inter-actuator coupling. 
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V. Conclusions 


This chapter presents a summary of the research completed as part of this Master’s 
work. Challenges overcome, key results, and future work are given. The research and 
observations made here build upon the previous work of Pellizzari and Schmidt, where this 
effort was aimed at transitioning their research towards a realistic implementation. In doing 
so, many new discoveries were made in the hopes of extending current AO capabilities in 
strong turbulence regimes. 

5.1 Challenges Overcome 

• Developed realistic AO system model. This was the crucial next step from the 
previous research. A realistic SRI model was developed, incorporating frame-by- 
frame fiber coupling effects, a first of its kind for this research. System radiometry 
was modeled, allowing for meaningful shot noise, senor read noise, and quantization 
noise effects. Care was given to correctly model the independent tilt and higher- 
order control systems, including the complete higher-order data processing pipeline. 
In doing so, data-flow optimizations were made specific to PCO-based wavefront 
reconstruction and AO control. 

• Complex, multi-dimensional parameter space. This research analyzed effects 
of Rytov number, Fried’s parameter, Greenwood frequency, and SNR for each 
wavefront reconstructor tested. This five dimensional parameter space created 
hundreds of gigabytes of simulation results which had to be efficiently tracked and 
managed. Running the simulations to explore the complete parameter space took 
weeks, even utilizing several high performance computers. Given the simulation 
requirements, significant development effort was given to optimize code. This 
framework is now available for future research. 
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• Bring AFIT AO hardware to an operational state. As this research began, AFIT’s 
AO system was plagued with a variety of issues that led to unreliable performance. 
Working in conjunction with the primary contractor, the author assisted in finding, 
diagnosing, and correcting a number of software and hardware issues. The author 
learned Lab VIEW in order to become proficient with the hardware interfacing and 
control code. Since no manuals were available, the author learned and wrote quick 
guides for common procedures. A full-system optical alignment procedure was 
developed, reducing the time required from 12 hours to only two. Unfortunately, 
at the final stages of this research a circuit board failure rendered the DM inoperable 
before any meaningful data could be taken. Initial AO system characterization results 
were presented at the 2011 SPIE Defense, Sensing, and Security conference. 

• Developed an ATS configuration change to allow high-Rytov simulations. In the 
configuration from the contractor, the ATS could simulate high-Rytov number 
conditions but only with prohibitively small r 0 ’s. In the simulations performed the 
ratio of DM actuator spacing to r 0 was kept, at most, to one. However, in the original 
configuration, the ATS could not maintain that ratio at cr 2 = 1. The author developed 
an optical modification that would allow much greater emulation capabilities in the 
r 0 : (t 2 x parameter space of interest to this research. 

5.2 Key Results 

• Demonstrated greatly improved strong turbulence AO performance in terms of both 
mean Strehl and Strehl variance. The AO performance gains using a PCO wavefront 
reconstructor over the traditional least squares reconstructor in terms of mean field 
estimated Strehl ratio and Strehl ratio normalized standard deviation are summarized 
in Table 5.1. Mean Strehl is a metric of absolute AO correction performance (higher 
is better), and the Strehl standard deviation is a metric of correction stability (lower 
is better). These results were presented at the 2012 IEEE Aerospace conference. 
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Table 5.1: Aggregate AO performance gains of all PCO-based wavefront reconstructors 
over the traditional least squares reconstructor. 


cr 2 x 

Mean Strehl 

Normalized Strehl Standard Deviation 

0.04 

0.2% 

-2.4% 

0.5 

19.1% 

-50.6% 

1 

126.4% 

-57.6% 


• Further reconstructor comparison under more realistic conditions. It turns out 
LSPV+1 is hard to beat in terms of performance gain relative computations required 
in most cases. However, as the Rytov number increases above roughly 0.5, a broader 
search becomes important as the PDF of h opt was shown to broaden. The simplicity 
of LSPV+4 with very similar performance to the more advanced reconstruction 
algorithms proposed by Pellizzari make it a good choice in high-Rytov conditions. 

• Identified cause of fades from PCO-based AO and demonstrated a potential solution. 
When in closed-loop operation in strong turbulence, PCO algorithms can cause 
short-term fades as h opt varies. These drops in Strehl are not caused by the higher- 
order control law as much as different tilt components associated with each PCO 
rotational realization. The tilt is removed from the higher-order corrections and 
instantly dumped on the tracking system; causing the fades. One solution was found 
by allowing tilt in the higher-order corrections with conservative controller gains. 
This prevented the fades with a slight tradeoff in mean Strehl ratio. Unfortunately, a 
realistic system with slight misalignments would require a zero-DC-gain control law 
to implement a similar solution. Additionally, this puts higher stroke requirements 
on the DM when tilt is being included. With this new understanding future work in 
control design could help mitigate these fades based on the application requirements. 
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Demonstrated and developed theory on improved PCO performance in weak 
turbulence, low SNR regimes. In weak turbulence, where no branch points are 
present, a consistent performance improvement of the PCO reconstructors over 
traditional least squares was found. Although the difference is small, this was 
completely unexpected based on previous concerns of noise corrupting the PCO 
process. A new theory proposed by the author shows in high-noise conditions, 
least squares reconstructors can introduce phase anomalies to the AO corrections. 
However, the PCO process removes these phase discrepancies introduced by the least 
squares reconstruction. These findings will be presented at the 2012 SPIE Optics and 
Photonics conference. 

Developed a set of theories regarding phase anomalies in high wavefront gradient 
conditions. When localized wavefront phases being reconstructed from an SRI are 
greater than n radians between adjacent subapertures, the least squares reconstructor 
will alias, leading to a phase anomaly in the reconstructed field. This can be mitigated 
using a higher density of SRI subapertures, at the cost of SNR. Alternatively, initial 
results show the PCO corrects the phase anomalies introduced by the least squares 
reconstruction. 

Formulated an optimized data flow for an AO system using PCO. Started by 
diagnosing coupling between the high-order and tracking control systems when using 
a PCO-based wavefront reconstructor. Since PCO reconstruction is not a single 
matrix operation, decoupling tilt from the higher-order commands required a new 
data flow process. The resulting process, illustrated in Fig. 3.4, minimizes matrix 
multiplies while maximizing AO effectiveness. 

Conducted first beam projection simulations using PCO wavefront reconstruction 
techniques. Conducted at the end of this research, this establishes a foundation 
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for future efforts to investigate these novel reconstruction methods with respect 
to beam projection applications. Although the initial simulation results identified 
needed changes, initial observations indicate promise over traditional least squares 
reconstruction. The identified flaws have since been fixed and follow-up simulation 
results will be presented at the 2012 DEPS Beam Control Conference. 

• Extensible AFIT AO system model for future research. The model is configurable 
to varying levels of realism and AO capabilities. Allows isolation of the research 
problem at hand without the difficulties of hardware work and the assumptions of 
analytical work. This includes the first modeling of realistic fiber coupling to our 
specific SRI optical setup - crucial to realistic noise modeling. 

5.3 Future Work 

• Hardware verification of simulation results. Unfortunately, the DM amplification 
circuit board failed prior to beginning hardware testing and verification of the 
simulation work presented here. This will be an essential step forward in verifying 
all theories formulated thus far regarding PCO-based wavefront reconstruction and 
control. 

• Hardware verifications using real atmospheric turbulence. Following hardware 
verification in the lab, the next step would be to test out the system performance 
using real atmospheric turbulence. As the realism is increased, new challenges are 
often discovered and additional validity is given to the theory. 

• Additional verification of phase anomaly theories. A set of initial simplified tests 
were made to form the theories presented on phase anomalies. More rigorous testing 
to relate the theories to more realistic scenarios is needed. 

• Further investigation into a practical tilt coupling between the control systems to 
reduce fades. A more in-depth consideration of the tilt coupling effects between 
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the tracking system and a higher-order system implementing a PCO wavefront 
reconstructor is needed. A zero-DC-gain controller is required to account for realistic 
component misalignments, which was not implemented here. Another alternative 
could be a feed-forward tilt-offloading between the two control systems as h is varied. 

Investigate adaptive controlling techniques with a PCO-based reconstructor. Promis¬ 
ing improvements in AO performance by Gibson [26] and Poyneer [32] using adap¬ 
tive AO control have been applied towards traditional weak turbulence techniques. 
Investigation as to whether these algorithms would work with the PCO non-unique 
rotational phase component would be needed. 

Effects of SRI performance with respect to platform vibrations. Moving the SRI out 
of a lab setup could reduce the fringe stability and thus the fidelity of the initial SRI 
wrapped measurements. This could be crudely simulated by disturbing the lab bench 
without requiring transportation of the system outside the lab. 

Further previous work on the hybrid SRI-Shack-Hartmann WFS using the simulation 
model created here. The simulation model developed as part of this research would 
be an excellent and easily extended platform on which to further previous work 
in designing an effective hybrid WFS. This may involve dynamic power splitting, 
comparative wavefront reconstruction, and more complex control techniques. The 
range of investigations and potential improvements are extensive, and this would be 
the first detailed test of a realistic implementation without the issues involved in full 
hardware integration. 

Further thought into IWCL and CL metrics in the PCO algorithms that operate 
over the residual phase error, as opposed to the fully aberrated field. In steady- 
state closed loop operation, the wavefront reconstructors do not observe the fully 
aberrated beacon phase. Rotational phase effects, including branch cuts, will be of 
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reduced visibility - often to the point where the branch cuts are below the threshold 
for detecting them at all. Only as the turbulence evolves in time are features visible, 
yet still not completely visible. These metrics do improve PCO performance with 
respect to the aliasing phase anomalies, but the author questions whether this is the 
best way to evaluate the “goodness” of a particular reconstruction when operating in 
closed loop. 
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Appendix A: Detailed Simulation Results 


The following plots show the compiled simulation results for each reconstructor over 
all conditions tested. Each grid-point represents a simulated one second AO run. Note that 
the simulations with higher Greenwood frequencies will run over many more Greenwood 
periods than the lower Greenwood frequency cases. Unfortunately, it would have been 
computationally prohibitive to run all of the low-Greenwood cases long enough to hold the 
number of Greenwood periods constant. While this is a definite compromise, I think the 
usefulness of the results is still strong. 

The left column is the mean field-estimated Strehl, calculated from after the higher- 
order control system is enabled through the end of each simulation run. The right column 
shows normalized standard deviation of the prior field-estimated Strehl ratio. This second 
statistic provides an idea for the system stability, where lower is better. The color axes have 
been fixed for each section below for easier comparison. 
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SNR 15 30 60 1 25 2 5 5 SNR 

Figure A.l: Least-Squares Reconstructor. Weak Turbulence: cr 2 = 0.04, d/r 0 = 1 


Mean Strehl: Normalized Standard Deviation 
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Figure A.2: LSPV+1 Reconstructor. Weak Turbulence: cr 2 = 0.04, d/r 0 = 1 
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Figure A.3: LSPV+4 Reconstructor. Weak Turbulence: cr£ = 0.04, d/r 0 = 1 


Mean Strehl: Normalized Standard Deviation 



Figure A.4: Histogram Reconstructor. Weak Turbulence: cr£ = 0.04, d/r 0 = 1 
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Figure A.5: Iterative Probing Reconstructor. Weak Turbulence: cr 2 = 0.04, <i/r 0 = 1 
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Figure A.6: Golden-Ratio Reconstructor. Weak Turbulence: cr 2 = 0.04, d/r 0 = 1 
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Figure A.7: Least-Squares Reconstructor. Medium Turbulence: cr 2 x = 0.5, d/r 0 = 1 
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Figure A.8: LSPV+1 Reconstructor. Medium Turbulence: cr 2 = 0.5, d/r 0 = 1 
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Figure A. 10: Histogram Reconstructor. Medium Turbulence: a 2 X = 0.5, d/r 0 = 1 
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Figure A. 12: Golden-Ratio Reconstructor. Medium Turbulence: cr 2 = 0.5, d/r 0 = 1 
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Figure A. 17: Iterative Probing Reconstructor. Strong Turbulence: cr 2 = 1, d/r 0 = 1 
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Appendix B: PCO Rotational Phase Tilt Propagation 


This appendix presents a simple simulation tracing a rotational wavefront through the 
higher-order data path. It is related to how the tracking system perceives the wavefront tilt 
both before and after the higher-order corrections are applied. 

The simulation starts with the rotational phase shown in Fig. B.l, with a uniform 
irradiance across the entire field. For this example, the wavefront sensor is idealized, 
giving no sensing errors or noise effects. The least squares reconstruction is also shown 
and is mostly flat as discussed in Sec. 2.3. Note that to the tracker, the rotational phase 
does not have any tilt. Figure B.2 shows the tracker spot for the uncorrected input field. It 
has a doughnut-like shape, but it is centered on the sensor array. 

The four PCO rotational reconstructions are plotted in Fig. B.3. In each rotational 
realization, a tilt component is defined by the location of the branch cut. It is important 
to distinguish this tilt component from the tilt component viewable by the tracker. The tilt 
component with respect to the higher-order system is defined by the tilt-removal operator 
from Sec. 3.1.1. To this operator, each rotational realization has a tilt component as shown 
in Fig. B.4. The tilt-removal operator is applied to all four rotational realizations, and the 
resulting outputs are applied to a DM. Sending the original rotational field onto each DM 
realization, and then looking again at what the tracker would see exposes an apparent shift 
in each of the tracker spots. Figure B.5 plots each tracker output. Note that the spots are 
no longer doughnut-like, so the higher-order phase aberrations were correctly compensated. 
However, the tilt components removed from each rotational realization, not only are visibly 
different, but they clearly propagate to the tracking system. Unfortunately, it does take a 
finite amount of time for the tracking system to react to these changes. That is the cause of 
the fades when h opt changes in a high Rytov environment. 
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Figure B.l: A hypothetical rotational phase (left), and its corresponding least squares 
reconstruction (right). 
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Figure B.2: The tracker image of the input field, from Fig. B.l. 
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Figure B.3: The four PCO rotational phase realizations from an LSPV+4 reconstructor 
over the rotational field from Fig. B. 1. 
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Residual Tilt, h = 0 


Residual Tilt, h = 7r/2 



Figure B.4: The four corresponding tilt components from each rotational phase realization 
of Fig. B.3. 
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Figure B.5: The four resulting tracker measurements after the higher-order corrections 
from each rotational phase realization of Fig. B.3. 
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